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Chapter  1 .  INTRODUCTION 


Evaluation  of  a  reconstruction  algorithm  should  be  done  using  a  sample  set  that  is  large  enough 
to  provide  us  with  a  statistically  significant  result  (Herman  and  Yeung,  1989).  In  order  to  carry  out 
an  evaluation,  one  possibility  is  to  use  a  set  of  computer  simulated  phantoms,  that  takes  into  account 
some  parameter  variabilities. 

This  technical  report  describes  in  detail  programs  that  generate  a  set  of  3D  phantoms  and  projection 
data,  reconstruct,  evaluate  and  then  compare.  The  main  characteristics  are: 

1.  Phantom  and  projection  data  generation: 

a.  Phantoms  with  many  (69)  ellipsoid  features,  ranging  from  small  (4  mm)  to  large  (40  mm) 
sized  features 

b.  Phantoms  are  random  samples  from  a  statistically  described  ensemble  of  3D  images  resembling 
those  to  which  PET  would  be  applied  in  a  medical  situation  (features  with  random  size, 
orientation  and  activity) 

c.  Features  are  inside  spheres  that  provide  background  value  for  some  important  clinical  tasks 
such  as  detectability 

d.  Types  of  features:  hot,  cold  and  normal  spots 

e.  Emulation  of  3D  PET  scanner  for  projection  data  generation,  with  detector  field-of-view  (FOV) 
blurring  and  a  realistic  3D  PET  noise  model 

2.  Reconstruction  algorithms: 

a.  Algebraic  Reconstruction  Technique  using  blob  (ARTblob)  as  basis  function  (Matej  and  Lewitt, 
1992) 

b.  ART  using  voxels  (ARTvox) 

c.  EM-ML  using  blobs  (EMblob) 

d.  EM-ML  using  voxels  (EMvox) 

3.  Evaluation  for  following  tasks: 

a.  Training  figure-of-merit  (FOM) 

b.  Structural  accuracy 

c.  Hot  spot  detectability 

d.  Cold  spot  detectability 

4.  Statistical  comparison  using  paired  t-test 

Justifications  for  using  some  models  (noise,  blurring,  contrast,  ...)  are  described  in  the  paper 
(Furuie  et  al,  1994)  and  an  application  evaluating  some  reconstruction  methods  is  reported  in  (Matej 
et  al,  1994). 

This  technical  report  also  describes  the  usage  and  assumptions  of  program  “genPhantom3d”,  that 
is  in  public  domain.  All  sources  are  provided,  basically  to  allow  users  to  compile  and  run  in  different 
computer  environments. 

These  programs  were  developed  using  C  language  (K&R),  under  a  UNIX  operating  system.  They 
have  been  tested  on  SUN  SPARCstations  (SunOS  Release  4,1.3)  and  Silicon  Graphics  machines  (IRIX 
Release  4.0.5  System  V),  however  we  disclaim  any  responsibility  for  possible  errors. 
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Chapter  2.  HOW  TO  INSTALL  AND  RUN 


2.1.  Requirements 

Disk  space  needed  for  sources  and  executable  programs  (static)  is  around  3  MB. 

For  the  default  configuration,  the  program  “genPhantom3d”  generates  at  least  3  files: 

1.  a  Phantom  Description  file  with  extension  .pde  (5,260  bytes) 

2.  a  phantom  file  with  extension  .img  (2,131,968  bytes) 

3.  a  projection  data  file  with  extension  .sen  (  33,939,456  bytes) 

However  a  file  that  contains  attenuation  factor  for  all  projection  rays  (phaLengths.AttenF.scn)  is 
also  required  in  order  to  implement  the  noise  model.  It  can  be  generated  once  simultaneously  with  the 
first  set  of  phantoms  (see  running  section).  Thus  the  minimum  amount  of  disk  space  is  around  70  MB. 


2.2.  Loading  the  installation  files 

The  tape  (1/4  inch,  150MB)  contains  two  tarfiles:  installation  files  and  3  samples  of  data  set. 
Procedure: 

1.  Create  a  directory  (e.g.  EVAL3DPET)  to  be  the  “evaluation”  directory  which  will  be  accessible 
to  users. 

2.  Change  directory  to  ’EVAL3DPET’  and  load  installation  files  from  tape  (tarfile),  using,  for 
example: 

tar  xvf  /dev/rstO  . 

This  procedure  will  create  and  load  the  following  directories: 

a)  sources 

b)  sources/gen 

gen.mak 

genPhantom3d.c 

MIPGgenPhan.c 

MIPGgenProjs.c 

MIPGphaFunc.c 

c)  sources/rec 

rec.mak 

ARTbloPet.c 

ARTvoxPet.c 

em_blpet.c 

em_vopeLc 

blob_fns.c 

basisfns.c 

besselfns.c 

voxeI_fns.c 
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d)  sources/bib 

bib.mak 

MIPGmath.c 

MIPGstring.c 

MIPGutiLc 

imagiol.l.c 

e)  sources/eval 

eval.mak 

eval3d.c 

complmgs.c 

student.c 

f)  sources/util 

utiLmak 

createPDE.c 

readPDE.c 

g)  include 

genPhantom3d.inc 

MIPGdefs.h 

MIPGgenPhan.h 

MIPGphaFunc.h 

MIPGgenProj  s.h 

MIPGstring.h 

MIPGmath.h 

MIPGutiLh 

blob_fns.h 

voxel_fns.h 

imagio.h 

pet_hdr.h 

h)  bin 

buildall.bat 

recon.bat 

eval.bat 

student.bat 

The  second  tarfile  contains  3  sets  of  phantoms  (phaMIPGOOO.*,  phaMIPGOOL*  and 
phaMIPGOlO.*,  where  *=.pde,  Jmg,  ^n)  and  attenuation  factor  file  (phaLengths.AttenF.scn). 
It  requires  around  145  MB.  It  will  create  directory  ’phantoms’. 
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2.3.  Compiling 

Once  the  installation  files  have  been  loaded,  there  is  a  script  that  allows  compilation  and  linking 
of  all  programs.  This  procedure  has  been  tested  on  the  SUN  and  Silicon  Graphics  workstations  and 
should  work  with  most  Unix  systems. 

Procedure: 

1.  Set  environment  variable  ’dir 3d’  to  the  “evaluation”  directory 
(Eg.:  setenv  dir3d  /usr/EVAL3DPET) 

2.  Change  directory  to  bin. 

3.  Run  the  script  buildall.bat,  for  example  with  parameter  ’spare’: 

buildalLbat  spare 

It  will  create  (if  not  already  there)  directory  spare  and  create  all  executables  there. 

2.4.  Running 

Line  commands: 

1.  Phantom/projection  generation: 

genPhantom3d  <prefix>  <start#>  <nuinPhantonis>  <totalCounts>  <FLAGatten> 

a.  prefix:  prefix  of  filenames  to  be  created.  Eg.:  phant 

b.  start#:  first  number  to  be  appended  to  prefix.  Eg.:  1(X) 

c.  numPhantoms:  number  of  phantoms  to  be  created.  Eg.:  5  (files  phantlOO.*  to  phantl04.* 
will  be  created) 

d.  totalCounts:  total  detected  counts  considering  all  detectors.  Eg.:  80000000 

e.  FLAGatten:  if  1,  it  will  also  create  file  with  attenuation  correction  factor 

(phaLengths.  AttenF.  sen) 

typical  CPU  (SUN  spare  10)  time  for  one  phantom/projection  data  generation  (complete,  with 
blurring  and  noise):  8  h. 

If  you  want  to  create  phantom  or  projection  data  for  different  reconstruction  volume,  phantom 
geometry  or  scanner  geometry  you  should  modify  the  file  “  genPhantom3d.inc”  (see  Appendix  C) 
prior  to  compilation.  This  file  contains  global  definitions  for  program  “genPhantom3d”. 

2.  Reconstruction  programs: 

For  all  following  reconstruction  programs,  it  is  assumed  that  the  input  projection  data  are  without 
attenuation.  Projection  file  (.sen)  contains  data  that  were  already  corrected  for  attenuation  (PET 
data  model),  or  contains  estimates  of  line  integrals  (CT  data  model).  (Time  information  is  for  1 
iteration,  using  SPARC  10,  see  Matej  et  al,  1994.) 


a.  ARTvox: 

Usage:  artvox  [Options  and  Parameters] 
Options: 

-h  for  this  message 
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Parameter  Values:  (note  space  before  value) 

-b  (Batch  mode  -  program  exits  if  wrong  parameters) 

-i  namepr  (name  of  input  projection  file) 

-f  name  (name  of  output^nput  data  &  image,  extensions  are  added) 

-s  startflag  (Start  from  :  l:data  on  disk;  2:const(avg)  cylinder,  3:zero) 

-n  noit  (number  of  iterations) 

-1  lambda  (relaxation  factor  (0.0  -  2.0)) 

reconstruction  time  (1  iteration):  2h  llmin 

b.  ARTblob: 

Usage;  artblob  [Options  and  Parameters] 

Options: 

-h  for  this  message 

Parameter  Values:  (note  space  before  value) 

-b  (Batch  mode  -  program  exits  if  wrong  parameters) 

-i  namepr  (name  of  input  projection  file) 

-f  name  (name  of  output/input  data  &  image,  extensions  are  added) 

-s  startflag  (Start  from  :  l:data  on  disk;  2:const(avg)  cylinder,  3:zero) 

-n  noit  (number  of  iterations) 

-1  lambda  (relaxation  factor  (0.0  -  2.0)) 

-r  blrad  (blob  radius) 

-m  blord  (order  of  Bessel  function  in  blob) 

-a  blalpha  (blob  parameter  alpha) 

(Typical/Standard  blob:  -r  2.0  -m  2  -a  10.4) 

reconstruction  time  (1  iteration):  10hl7min 

c.  EMvoxel: 

Usage:  emvox  [Options  and  Parameters] 

Options: 

-h  for  this  message 

Parameter  Values:  (note  space  before  value) 

-b  (Batch  mode  -  program  exits  if  wrong  parameters) 

-i  namepr  (name  of  input  projection  file) 

-f  name  (name  of  output/input  data  &  image,  extensions  are  added) 

-s  startflag  (Start  from  :  l:data  on  disk;  2:constant  cylinder) 

-n  noit  (number  of  iterations) 

-t  nameattpr 

(filename  -  of  attenuation  coefficients  on  proj. lines  (program  attenuates  projections  and 
considers  attenuation  also  in  EM-ML  algorithm), 

0  -  if  CT  or  no  attenuation  of  PET  data  considered.) 

[  -c  comame  (base  name  of  corr.  matrix  file  -  if  other  than  standard)] 
reconstruction  time  (1  iteration):  3hl0min 
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d.  EMblob 

Usage:  emblob  [Options  and  Parameters] 

Options: 

-h  for  this  message 

Parameter  Values:  (note  space  before  value) 

-b  (Batch  mode  -  program  exits  if  wrong  parameters) 

-i  namepr  (name  of  input  projection  file) 

-f  name  (name  of  output/input  data  &  image,  extensions  are  added) 

-s  startflag  (Start  from  :  l:data  on  disk;  2:constant  cylinder) 

-n  noit  (number  of  iterations) 

-r  blrad  (blob  radius) 

-m  blord  (order  of  Bessel  function  in  blob) 

-a  blalpha  (blob  parameter  alpha) 

(Typical/Standard  blob:  -r  2.0  -m  2  -a  10.4) 

-t  nameattpr 

(filename  -  of  attenuation  coefficients  on  proj. lines  (program  attenuates  projections  and 
considers  attenuation  also  in  EM-ML  algorithm), 

0  -  if  CT  or  no  attenuation  of  PET  data  considered.) 

[  -c  comame  (base  name  of  corr.  matrix  file  -  if  other  than  standard)] 
reconstruction  time  (1  iteration):  9h  lOmin 


For  programs  ’artvox’,  ’artblob’,  ’emvox’,  ’emblob’,  one  can  see  the  running  status  by  printing 
file  ’name’.stat.  The  file  ’name’. hist  gives  information  about  parameters  used  for  reconstruction. 
Each  of  these  4  programs  can  be  used  sequentially.  For  example,  in  order  to  use  ’emvox’  and 
store  results  at  10th  and  30th  iteration,  with  starting  image  as  a  constant  cylinder  (-s  2),  one  could 
use  the  following  script: 

#  Reconstruct  first  10  iterations  and  evaluate 

emvox  -b  -i  phantom. sen  -f  phantomEmvox  -s  2  -n  10  -t  0 
eval3d  phantom.img  phantomEmvox.img  phantomEmvox.eval  4 

#  Reconstruct  next  20  iterations  and  evaluate 

emvox  -b  -i  phantom.scn  -f  phantomEmvox  -s  1  -n  20  -t  0 
eval3d  phantom.img  phantomEmvox.img  phantomEmvox.eval  4 
cat  phantomEmvox.hist  »  phantomEmvox.eval 

All  these  4  programs  assume  that  there  is  a  .pde  file  associated  to  projection  data  file.  In  the 
example  above,  it  is  assumed  that  ’phantom.pde’  exists.  It  is  possible  to  create  a  ’.pde’  file  using 
program  createPDE  (see  item  4). 

3.  Evaluation: 
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a.  eval3d  <ideal  image  filename>  <reconstructed  image  filename>  <evaluation  result  filename> 
<DETEC'I^  [FLAG] 

where: 

DETECT:  detectability  index 
0:  equation  (4)  [see  section  3.3] 

4:  equation  (5)  [see  section  3.3] 

FLAG  (optional): 

0:  default  (prints  only  numerical  values) 

1:  prints  ROC  global  curve 
2:  prints  ROC  global  &  partial  curves 

eval3d  utilizes  .pde  file  of  <ideal  image  filename>  to  get  information  for  evaluation  (dimen¬ 
sionalities,  position,  size  and  activity  of  features,...). 

This  program  appends  (creates  if  not  exist)  results  of  evaluation  to  the  file  <evaluation  result 
filename>.  In  order  to  facilitate  further  analysis  (such  as  statistical  comparison)  it  also  cre¬ 
ates  a  file  <evaluation  result  filename>.tbl  that  contains  (cumulatively)  a  table  of  calculated 
figure-of-merit  (FOMs).  For  each  call  to  evaI3d  it  stores  in  the  .tbl  a  row  with: 
<reconstructed  image  filename>  cTraining  FOM>  <SA_gl>  <HSD_gl>  <CSD_gl>  <SA_S> 
<SA_M>  <SA_L>  <HSD_S>  <HSD_M>  <HSD_L>  <CSD_S>  <CSD_M>  <CSD_L> 
<HSD_avg>  <CSD_avg> 
where  Training  FOM:equation  7, 

SA:  structural  accuracy  (equation  3), 

HSD:  hot  spot  detectability  using  ROC  area, 

CSD:  cold  spot  detectability  using  ROC  area. 

The  suffix  means: 

_gl  :  global  value,  considering  features  of  all  sizes, 

_avg:  averaged  ROC  area,  calculating  separately  ROC  area  for  small  (_S),  medium  (_M) 
and  large  CL)  sizes, 

_S:  small  features  (surrounding  sphere  radius  [4.0 , 6.4}  mm), 

_M:  medium  size  features  (surrounding  sphere  radius  (8.0  ,  10.0}  mm), 

_L:  large  size  features)  (surrounding  sphere  radius  (16.0 , 20.0}  mm). 

There  are  high  similarities  among  some  of  FOMs  (see  appendix  G)  for  instance  between 
HSD_gl  and  HSD_avg  using  equation  4  or  between  HSD_avg  using  equation  4  and  5. 

b.  complmgs  <reference.img>  <imageToBeCompared.img> 

It  calculates  minimum,  maximum,  average  and  total  counts  for  <reference.img>,  «dmageTo- 
BeCompared>  and  (<reference.img>  —  <imageToBeCompared>). 

c.  student  <fileTbll>  <fileTbl2>  <numbRows>  <numCols>  <column[0..]>  [<FLAG>] 

It  compares  (student  statistical  test)  columns  (<column>)  of  two  tables  (or  same  table)  in  files 
<fileTbll>  and  <fileTbl2>,  each  with  <numbRows>  rows  and  <numCols>  columns.  [FLAG] 
is  optional  and  if  equal  to  1,  it  prints  details  of  calculation.  See  Appendix  F  for  examples. 

4.  Others: 

a.  createPDE 
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It  is  a  simple  program  to  create  a  .pde  file.  Sometimes  we  have  a  projection  data  (such  as 
Hoffman  phantom)  in  Petview  format  and  we  want  to  reconstruct  it  using  above  mentioned 
iterative  methods.  In  this  case  we  have  to  create  a  .pde  file. 

b.  readPDE 

It  reads  and  displays  the  content  of  a  .pde  file 

c.  recon.bat,  eval.bat,  student.bat 

Samples  of  scripts  (see  Appendix  F)  for  reconstruction,  evaluation  and  statistical  comparison 
in  batch. 
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Chapters.  FRAMEWORK 


In  this  chapter  we  describe  the  phantom  model,  projection  data  generation,  some  figures  of  merit 
for  evaluation  and  testing.  For  a  complete  description  with  justifications,  see  (Furuie  et  al,  1994). 

3.1.  Phantom  generation 

In  order  to  evaluate  fully  3-D  PET  reconstruction  algorithms,  we  are  using  an  ensemble  of 
phantoms,  which  are  roughly  related  to  the  task  of  imaging  features  of  relatively  high  or  low  activity 
(as  compared  to  their  background)  in  the  head.  The  “head”  is  an  ellipsoid  (with  radii  Rx  =  128.0 
mm,  Ry  =  96.0  mm,  Rz  =  147.8  mm,  and  center  at  z  =  64  mm,  see  Fig.  1),  all  of  which  is  assumed 
to  contain  attenuating  material,  but  only  a  part  of  which  (the  part  which  lies  in  the  “reconstruction 
region”  indicated  in  Fig.  1)  is  assumed  to  have  activity  in  it.  The  activity  in  that  part  is  assumed 
to  be  uniform  (and,  since  we  are  free  to  choose  our  units,  we  consider  this  uniform  activity  to  have 
value  1)  except  for  69  ellipsoid-shaped  features.  Although  the  features  themselves  change  in  size  and 
shape  from  phantom  to  phantom  (as  is  explained  below),  each  of  the  features  is  located  inside  one 
of  69  spheres,  which  themselves  are  the  same  in  all  the  phantoms.  The  geometrical  distribution  of 
the  spheres  is  indicated  in  Fig.  1  and  their  sizes  are  shown  in  Table  I.  Precise  location  and  size  of 
spheres  are  given  in  Appendix  D.  As  indicated  on  the  right  side  of  Fig.  1,  the  spheres  are  located 
in  three  layers  (A,  B,  and  C).  Each  layer  contains  11  larger  spheres,  whose  centers  are  distributed 
according  to  the  schematic: 
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In  the  middle  of  the  segments  that  connect  the  center  of  sphere  ‘04’  to  that  of  ‘01’  or  of  ‘09’  and  that 
connect  the  center  of  sphere  ‘08’  to  that  of  ‘03’  or  of  ‘11’,  there  are  smaller  spheres  (see  left  side  of 
Fig.  1).  Immediately  above  and  below  (in  the  z-direction)  each  of  these  smaller  spheres  there  are  two 
other  spheres  of  the  same  radius  and  so  we  have  altogether  12  smaller  spheres  per  layer.  The  centers 
of  the  layers  C,  B,  and  A  are  at  heights  -50  mm,  -16  mm,  and  34  mm,  respectively. 

The  features  are  generated  using  the  following  random  process: 

•  The  shape  of  a  feature  is  an  ellipsoid  whose  centre  is  at  the  center  of  one  of  the  spheres  discussed  in 
the  previous  paragraph.  Each  of  the  three  radii  of  the  ellipsoid  is  independently  and  randomly  selected 

Table  I.  Radii  of  the  spheres  containing  features,  contrast  of  the  features  in  them  if  they  are 
“hot”  or  “cold,”  their  location,  and  the  total  number  of  spheres  of  each  type  in  a  phantom. 


Radius 

Contrast 

Layer 

Number 


4.0  mm 
75.0% 
C 
12 


6.4  mm 
29.3% 
B 
12 


8.0  mm 
18.8% 
A 
12 


10.0  mm 
12.0% 
C 
11 


16.0  mm 
4.7% 

B 

11 


20.0  mm 
3.0% 

A 

11 
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256  mm 


Reconstruction  region 


a> 
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Transverse  slice  (at  2=34  mm) 


Coronal  slice  (at  y=0  mm) 


Figure  1.  Illustration  of  the  shape  of  the  reconstruction  region  in  our  phantoms  and  of 
the  geometrical  distribution  of  the  spheres  which  contain  the  randomly-shaped  features. 


from  a  distribution  which  is  uniform  between  50%  to  100%  of  the  radius  of  the  enclosing  sphere. 
Orientations  of  the  ellipsoid  axes  (rotations  around  the  three  coordinate  axes)  are  also  independently 
and  randomly  selected  from  a  distribution  which  is  uniform  between  0  and  180  degrees.  The  part 
of  the  enclosing  sphere  which  is  outside  the  feature  is  referred  to  below  as  the  “background”  of  the 
feature.  (Note  that  since  we  have  features  of  varying  sizes,  our  methodology  can  be  adapted  to  allow 
us  to  evaluate  reconstructions  from  the  points  of  view  of  FOMs  specialized  to  features  in  a  certain 
size  range.) 

•  Each  feature  is  randomly  assigned  to  be  one  of  three  equiprobable  types:  NORMAL  (the  activity  in 
the  feature  is  the  same  as  in  its  background),  HOT  (the  feature  is  a  “hot  spot”),  or  COLD  (the  feature 
is  a  “cold  spot”).  The  activity  assigned  to  hot  spots  (respectively,  to  cold  spots)  is  higher  (respectively, 
lower)  than  the  activity  in  their  background  by  the  contrast  listed  in  Table  I.  This  contrast  varies  from 
3%  to  75%  (of  the  background  activity),  depending  on  the  size  of  the  enclosing  sphere.  Such  variation 
in  contrast  is  necessary  in  order  to  have  approximately  the  same  difficulty  in  detecting  features  with 
different  sizes  (see  Appendix  A).  The  actual  size  of  the  contrast  (only  3%  for  large  features  and 
75%  for  very  small  features)  is  much  less  than  what  one  would  like  to  have  for  reliable  detection  of 
features.  However,  this  is  essential  for  comparing  reconstruction  algorithms;  if  all  features  had  high 
contrast,  then  we  would  have  nearly  100%  detectability  associated  with  all  reasonable  methods  and  we 
would  not  be  able  to  determine  their  relative  merit  for  hot  spot  and  cold  spot  detection.  (We  do  not 
include  a  phantom  in  our  ensemble  unless  it  has  at  least  one  feature  of  each  type.)  For  the  evaluation 
computation,  we  consider  that  a  voxel  is  “in  a  feature,”  if  its  center  and  at  least  five  of  its  vertices 
are  in  the  feature  (i.e.,  they  lie  in  the  ellipsoid  which  is  the  support  of  the  feature).  We  consider  that 
a  voxel  is  “in  the  background  of  a  feature,”  if  its  center  and  at  least  five  of  its  vertices  are  in  the 
sphere  enclosing  the  feature  but  not  in  the  feature.  If  the  random  generation  of  a  phantom  resulted  in 
it  containing  a  feature  such  that  either  the  feature  itself  or  its  background  contained  fewer  than  five 
voxels,  then  we  generate  another  feature  and  check  again. 
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Figure  2.  Slices  through  one  of  the  randomly  generated  phantoms 
(first  row:  transverse;  second  row:  sagittal;  third  row:  coronal). 

For  comparisons  with  reconstructions,  each  phantom  is  digitized  as  an  array  with  64x128x128 
cube-shaped  voxels,  each  with  2mm  side  length.  The  value  assigned  to  a  voxel  in  a  phantom  is 
determined  as  the  average  of  the  values  of  activity  at  the  regularly-spaced  points  of  a  4x4x4  array 
within  the  voxel.  Slices  through  the  digitized  version  of  one  of  our  randomly  generated  phantoms 
are  shown  in  Fig.  2. 

3.2.  Projection  data  generation 

For  the  generation  of  projection  data,  we  simulate  an  actual  3-D  PET  scanner  with  a  cylindrical 
detector  of  420  mm  in  diameter  and  256  mm  in  axial  length;  see  Fig.  3.  (Such  a  scanner  is  currently 
under  construction  in  our  department.)  We  are  assuming  that  data  are  rebinned  into  sets  of  parallel 
raysums,  so  that,  for  each  “view,”  we  have  a  virtual  panel  of  90  x  128  measurements  corresponding  to 
“bins”  perpendicular  to  the  panel.  This  panel  “rotates”  180  degrees  around  the  reconstruction  volume 
producing  a  total  of  96  views  per  tilt  position.  We  generate  15  steps  of  panel  tilting  between  —26.5  to 
26.5  degrees  (step  =  3.75  degrees).  It  is  assumed  that  the  center  line  (axis  for  tilting)  of  the  panel  is  in 
the  central  plane  of  the  reconstruction  volume  (z=0)  and  contains  the  origin  of  coordinate  system.  Thus 
the  total  number  of  bins  is  15x96x90x128  (16.6  million).  The  cylindrical  nature  of  om  simulated 
detector  implies  that  raysums  are  not  available  for  some  of  the  bins  in  the  tilted  views.  In  calculating 
the  projection  data  for  those  bins  for  which  such  data  are  available,  the  cross-section  (2x2  mm)  of 
the  bin  is  sampled  by  a  square  arrangement  of  four  rays  (parallel  to  and  1  mm  apart  from  each  other). 
For  each  of  these  rays  the  integral  of  activity  is  calculated  analytically.  These  calculated  integrals  are 


420  mm 


Figure  3.  Illustration  of  scanner  geometry.  The  reconstruction  volume  is  centered  in  the  scanner. 

blurred  by  a  Gaussian  smoothing  filter  with  a  full-width  at  half-maximum  (FWHM)  approximately 
equal  to  4  mm,  to  match  the  combined  effects  of  the  uncertainty  in  the  exact  location  of  a  detected 
event  and  of  the  rebinning  from  actual  measurements  by  such  a  scanner  into  the  arrangement  of  panels 
that  we  use  in  the  reconstructions.  The  blurred  integrals  which  are  associated  with  the  four  rays  of  a 
bin  are  averaged  to  provide  us  with  the  “raysum”  of  activity  associated  with  that  bin. 

Because  of  attenuation,  the  measurement  (detected  counts)  associated  with  a  bin  will  be  substan¬ 
tially  less  than  its  raysum.  We  take  the  expected  value  of  the  detected  counts  for  a  particular  bin  to 
be  pik,  where  p  is  the  raysum  associated  with  the  bin  and  the  attenuation  factor  k  is  calculated  to  be 
expiLp),  where  L  is  the  analytically  calculated  length  of  intersection  (in  mm)  of  the  central  ray  of  the 
bin  in  question  with  the  phantom  and  p  is  the  linear  attenuation  coefficient  of  water  (0.0095  mm“^ 
for  511  keV).  For  our  simulations,  we  assumed  the  total  number  of  detected  counts  to  be  80  million, 
which  represents  nearly  56,000  per  projection  plane  and  averages  to  be  around  5  per  bin.  The  actual 
range  of  the  expected  number  of  detected  counts  per  bin  for  our  phantoms  turns  out  to  be  from  0  to  27. 
After  attenuation  correction,  the  total  number  of  counts  in  sinogram  (all  panels)  is  around  335  million. 

Our  calculation  of  the  noise  in  the  projection  data  takes  into  consideration  that  corrections  for 
attenuation,  scattering,  detector  normalization,  etc.  cause  a  discrepancy  between  the  noise  in  actual 
data  and  the  noise  that  would  be  provided  by  a  pure  Poisson  model  (see,  e.g.,  Rowe  and  Dai  1992). 
For  any  calculated  raysum  p,  we  replace  it  by  its  noisy  version,  which  is  calculated  as 

kx  Poisson -I- TV (0,5“)  ,  (1) 

where  Poisson(plk)  denotes  a  sample  from  a  Poisson  distributed  process  with  mean  pIk  (the  attenuated 
raysum)  and  ^^(O,  s^)  denotes  a  sample  from  a  zero-mean  normal  distribution  with  standard  deviation 
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equal  to  s  (whose  actual  value,  as  we  show  below,  depends  on  pik  in  a  well-defined  way).  This  second 
random  distribution  forms  part  of  our  noise  model  to  account  for  discrepancies  between  actual  physical 
measurements  and  our  model  of  them.  The  methodology  for  choosing  s  in  conjunction  with  the  noise 
model  implied  by  the  form  of  equation  (1)  was  determined  experimentally  using  cylinder-shaped 
physical  phantoms  (Furuie  et  al,  1994): 


Our  projection  data  generation  process  can  be  summarized  as  follows:  raysums  (p)  of  activities 
are  first  generated  using  analytical  integrations  of  the  geometrically  described  phantom  and  a  blurring 
process  and  this  is  then  followed  by  the  introduction  of  noise  according  to  equations  (1)  and  (2). 
Samples  from  one  of  the  resulting  projection  data  sets  are  shown  in  Fig.  4  for  tilts  at  0  and  at  -26.5 
degrees.  For  the  latter  case,  we  also  see  the  effects  of  not  having  measurements  for  a  large  number 
of  bins  due  to  the  cylindrical  nature  of  our  detector. 


Figure  4.  Samples  of  projection  data  for  views  with  tilt  of  0  degrees  (top)  and  with  tilt  of  — ^26.5 
degrees  (bottom).  Observe  that  data  are  not  available  for  quite  a  few  bins  in  the  latter  case 

3.3.  Figures  of  merit 

The  definition  of  a  figure  of  merit  (FOM)  should  be  related  to  the  task  of  solving  a  medical 
problem.  Prior  to  giving  some  examples,  we  discuss  some  notation  common  to  all  of  them.  In  this 
notation  we  assume  a  single  phantom  and  a  single  reconstruction  (by  some  method)  from  a  projection 
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data  set  of  that  phantom.  An  FOM  is  a  measure  of  how  good  that  particular  reconstruction  is  from 
some  mathematically  precise  point  of  view. 

Consider  a  phantom  which  contains  F  features.  (In  all  our  phantoms  F  is  69,  but  the  definitions 
of  the  FOMs  are  applicable  to  a  much  larger  class  of  phantoms.)  We  use  f  to  refer  to  the  fth  of  the 
F  features.  For  this  particular  feature,  we  use  Rf  to  denote  the  radius  of  the  surrounding  sphere,  fif 
to  denote  the  average  activity  of  the  voxels  within  the  feature  in  the  digitized  phantom,  nf  to  denote 
the  number  of  voxels  in  the  feature,  m/  and  vy  to  denote  the  average  and  the  variance  of  the  activity 
within  the  feature  in  the  reconstruction,  n/  to  denote  the  number  of  voxels  in  the  background  of  the 
featme,  and  my  and  uy  to  denote  the  average  and  the  variance  of  the  activity  within  the  background 
of  the  feature  in  the  reconstruction. 

We  define  the  FOM  structural  accuracy  of  a  reconstruction  as  the  value  of 

(3) 

where  we  think  of  the  normalizing  parameter  Rq  as  the  radius  of  a  reference  sphere.  Specifically, 
for  our  phantoms  we  select  Rq  to  be  9.0  mm,  which  is  an  intermediate  size  of  our  range  of  radii 
(see  Table  I).  Note  that  more  accurate  average  activity  values  in  the  features  in  the  reconstruction 
will  result  in  a  higher  value  of  structural  accuracy,  with  a  perfect  reconstruction  yielding  the  value  1 
for  this  FOM.  (For  our  phantoms  we  also  found,  see  Matej  et  al  1994,  that  practical  reconstruction 
algorithms  yield  structural  accuracy  values  in  excess  of  0.5  and  so  in  the  same  range  as  provided  by 
other  FOMs  that  we  discuss  below.)  The  sum  in  (3)  is  a  “weighted  1-norm”  for  the  average  activities 
within  features,  where  the  weights  account  for  the  different  sizes  of  the  features.  These  weights  were 
derived  from  the  following  considerations.  If  we  assume  that  the  variance  of  the  noise  in  a  single 
voxel  of  the  reconstruction  is  the  same  for  all  voxels,  then  the  variance  of  the  noise  in  my  is  inversely 
proportional  to  wy  and  so,  on  average,  to  R^.  Therefore,  under  the  circumstance  that  in  each  individual 
feature  the  standard  deviation  of  noise  is  the  same  for  all  voxels  in  the  feature,  our  chosen  weight 
factors  would  yield  a  weighted  1-norm  independent  of  size,  which  is  proportional  to  the  average  of  the 
standard  deviations  of  error  in  the  individual  features.  The  value  of  the  unweighted  1-norm  would  be 
influenced  mainly  by  the  noise  in  the  smaller  features;  by  raising  the  normalizing  factor  to  the  power 
3/2  we  equalize  the  expected  influence  of  features  of  various  sizes. 

For  the  next  two  FOMs  (Hot  Spot  Detectability  and  Cold  Spot  Detectability)  we  need  an 
intermediate  concept.  Wfith  every  feature  /  in  a  reconstruction  we  associate  a  detectability  index 
hf.  One  possibility  is  to  use  index 


which  is  based  on  one  of  the  standard  measures  of  separability  (between  the  reconstructed  voxel  values 
in  the  feature  and  the  reconstructed  voxel  values  in  the  background  of  the  feature)  in  statistical  pattern 
recognition  theory  (Fukunaga  1972)[see  Appendix  B  for  details]. 

Note  that  in  a  good  reconstruction,  we  would  expect  this  index  to  be  positive  if  the  feature  is  a 
hot  spot,  negative  if  the  feature  is  a  cold  spot,  and  approximately  zero  otherwise.  The  absolute  value 
of  the  index  goes  up  if  the  absolute  difference  in  the  means  of  voxel  values  in  the  feature  and  in  its 
background  gets  greater,  but  goes  down  if  the  variance  in  these  voxel  values  goes  up. 
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In  the  evaluation  of  the  practical  reconstructions  there  were  several  problems  with  the  detectability 
index  (4)  derived  on  the  basis  of  idealistic  theoretical  considerations.  First,  in  order  to  obtain 
comparable  detectability  for  the  features  of  different  sizes,  the  feature  contrast  has  to  be  changed 
more  than  suggested  by  the  formula  (23),  which  was  considered  in  derivation  of  index  (4).  This 
is  caused  by  reconstruction  problems,  like  border  effects,  which  are  dependent  on  the  feature  size. 
Secondly,  since  noise  in  the  reconstruction  is  correlated  on  the  local  neighborhood,  we  are  not  able  to 
accurately  estimate  the  local  variances  (v/,  for  small  features.  Another  possibility  is  to  approximate 
them  by  global  variance  for  given  reconstructed  image.  In  this  case,  the  variances  (v/,  will  have 
the  same  value  for  all  features  in  the  given  reconstruction.  Therefore,  their  value  will  not  influence 
the  ROC  analysis  result  and  they  can  be  omitted  from  the  detectability  index.  Thirdly,  in  the  practical 
reconstructions  using  several  reconstruction  algorithms,  it  has  been  found  that  the  size  dependent 
correction  factor,  in  detectability  index,  given  by  the  contrast  change  does  not  have  desirable  effect 
Practically,  it  has  been  found  that  the  best  results  (the  closest  approximation  to  the  independence  of 
expected  value  and  variance  of  hf  on  the  size  of  features)  were  obtained  using  the  hf  without  any 
correction  at  all.  Considering  second  and  third  points  above,  the  best  results  were  obtained  in  our 
practical  evaluations  using  following  detectability  index; 


(5) 


This  detectability  index  turned  out  to  be  less  sensitive  to  feature  size  for  the  given  ensemble  of 
phantoms.  Therefore,  equation  5  is  preferable  for  ROC  analysis  using  our  phantoms.  The  evaluation 
program  eval3d  provides  ROC  area  considering  both  detectability  indices  and  also  ROC  area  per  size 
(see  section  2.4.3). 

We  define  the  FOM  hot  spot  detectability  of  a  reconstruction  as  the  area  under  the  ROC  (Receiver 
Operating  Characteristic)  curve  based  on  the  index  hf  (Swets  and  Pickett  1982,  Green  and  Swets  1988, 
Herman  and  Yeung  1989).  To  be  exact,  we  do  the  following.  Suppose  that  there  are  M  hot  spots 
and  N  normal  features  in  the  phantom  under  consideration.  For  any  real  number  h,  let  c{h)  denote 
the  number  of  HOT  features  /  for  which  hf>h.  Let  be  the  NORMAL  features,  with  the 

notation  assigned  so  that  hj^  >  hf^_^  >  hf^^.  Then  hot  spot  detectability  is  defined  to  be 


MN  ^  ■ 

n=l 

It  is  easy  to  see  that  this  FOM  has  a  value  between  0  and  1  and  that  the  value  gets  near  the  high  end  if 
the  indices  hf  for  normals  tend  to  be  small  as  compared  to  those  for  hot  spots.  In  fact,  for  any  reasonable 
reconstruction  method,  we  expect  the  hot  spot  detectability  to  be  greater  than  0.5.  (A  method  which 
randomly  assigns  values  to  voxels  should  end  up  with  hot  spot  detectability  approximately  0.5.) 

Except  for  obvious  minor  changes,  cold  spot  detectability  of  a  reconstruction  is  defined  along  the 
same  lines. 

The  three  FOMs  defined  so  far  have  been  designed  to  reflect  various  medical  tasks:  structural 
accuracy  is  a  measure  of  how  well  we  can  estimate  the  total  uptake  in  features  and  the  other  two 
indicate  how  well  we  can  detect  hot  spots  and  cold  spots,  respectively.  We  complete  our  list  of  FOMs 
with  a  measure  which  is  less  relevant  from  the  medical  point  of  view,  but  which  we  have  found  useful 
for  training  algorithms  (see  the  next  section  and  Matej  et  al  1994).  For  this  reason  we  call  it  the 
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training  FOM.  It  is  defined  as 


where  ef  is  the  average  over  voxels  in  the  yth  feature  of  the  absolute  difference  between  voxel  values 
in  the  reconstruction  and  in  the  phantom. 

3.4.  Training  and  testing 

Many  reconstruction  algorithms  have  some  free  parameters  in  them  (such  as  relaxation  parameters 
and  the  number  of  iterations  in  iterative  methods  or  filter  characteristics  in  transform  methods,  Herman 
1980).  To  insure  that  a  general  reconstruction  method  which  has  some  free  parameters  is  not  unjustly 
discarded  as  a  result  of  a  bad  choice  of  these  free  parameters  in  an  evaluation  study,  it  is  appropriate 
that  before  testing  the  method  one  finds  those  values  of  the  free  parameters  for  which  it  performs 
optimally  (or  at  least  well)  according  to  some  relevant  FOM.  This  is  referred  to  as  training.  In  our 
environment  training  can  be  done  as  follows. 

We  randomly  generate  from  our  ensemble  a  few  phantoms  and  their  associated  projection  data; 
we  refer  to  this  as  the  training  set.  We  then  search  for  values  of  the  free  parameters  which  optimize 
the  training  FOM  for  the  training  set.  (In  practice,  we  limit  our  search  until  we  find  values  of  the  free 
parameters  which  are  such  that  analysis  of  changes  in  the  training  FOM  with  changing  free  parameters 
indicate  that  we  are  within  a  few  percent  of  the  optimal  value.)  The  frree  parameters  then  are  fixed 
at  these  values  for  use  in  the  testing  process  described  below.  It  is  reasonable  in  this  numerical 
observer  approach  to  use  a  different  FOM  for  training  and  for  testing.  This  is  because  the  performance 
of  algorithms  can  sometimes  be  fine-tuned  to  a  very  high  level  by  the  use  of  training  on  the  same 
FOM  as  is  used  for  testing  (Herman  and  Yeung  1989).  This  is  likely  to  be  unrealistic  for  a  clinical 
application,  where  we  would  not  have  the  luxury  of  being  able  to  use  a  different  most-appropriately 
optimized  reconstruction  algorithm  for  each  of  the  different  tasks. 

The  evaluation  should  be  performed  using  a  testing  set  of  phantoms  and  projection  data  which  is 
fix)m  the  same  ensemble  as,  but  is  statistically  independent  of,  the  training  set.  Each  reconstruction 
algorithm  that  is  being  compared  should  be  applied  to  all  the  projection  data  and  FOMs  should  be 
calculated  for  each  of  the  reconstructions. 

Once  an  FOM  is  calculated  for  each  reconstruction  produced  by  two  different  reconstruction 
algorithms,  we  assign  a  level  of  statistical  significance  to  rejecting  the  null-hypothesis  that  the  two 
reconstruction  methods  are  equally  good  (from  the  point  of  view  of  the  FOM),  in  favor  of  the  hypothesis 
that  the  one  with  the  higher  FOM  is  better,  as  follows.  Let  /3b  and  6b  be  the  FOMs  of  the  reconstructions 
from  the  6th  of  a  total  of  B  data  sets,  as  reconstructed  by  the  two  methods,  respectively.  Then, 
according  to  the  null-h5q)othesis,  /3b  —  6^  is  a  sample  of  a  zero-mean  random  variable.  It  follows, 
for  large  enough  B,  that 

B 

Yli^b-Sb)  (8) 

6=1 

is  a  sample  of  a  normally  distributed  zero-mean  random  variable  (Mould  1989).  The  variance  of  this 
random  variable  is  B  times  that  of  the  zero-mean  random  variable  of  which  ^b  —  h  are  samples  for 
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enough  B,  it  is  reasonable  to  assume  that  the  null-hypothesis  impUes 
that  (8)  is  a  sample  from  a  normally  distributed  random  variable  with  mean  zero  and  variance 

(^6  -  hf  • 

b=i 

We  can  thus  use  the  normal  distribution  to  calculate  significance  (Mould  1989). 

An  alternative  way  of  assigning  significance,  which  may  be  preferred  when  B  is  considered  not 

arge,  is  by  using  the  r-test  for  paired  data  with  the  value  of  t  calculated  according  to  the  following 
formulas  (Mould  1989): 

.  m 

*^7J7S’ 

Ei^b-Si) 

m  =  - ,  and  (10) 


S  ((/^i  -  ^b)  -  rnf  . 
tel 

The  calculated  t  (eq.  10)  was  compared  to  Student  table  for  DF  (degree  of  freedom)  equal  to  (5-7), 
for  a  statistical  level  of  significance  (5%,  0.5%  and  0.05%). 

3.5.  Implemented  Reconstruction  Algorithms 

The  Iterative  reconstruction  algorithms  may  be  perceived  as  attempts  to  solve,  in  some  sense,  the 
system  of  equations 


y  =  Ax, 

where  y  is  the  vector  of  projection  data,  A  is  a  matrix  whose  entries  are  the  contributions  of  single 
image  basis  elements  to  the  data  associated  with  particular  bins,  and  x  is  the  vector  of  the  unknown 
mage  values  (coefficients  in  the  decomposition  of  the  image  into  the  basis  elements).  Traditionally,  the 
discrete  model  of  an  iterative  reconstruction  method  assumes  voxel  basis  functions  (i.e.,  basis  function 
Whose  value  is  1  inside  a  specific  cube-shaped  voxel  and  is  0  everywhere  else).  An  alternative  approach 
IS  discussed  in  (Lewitt  1990,  Lewitt  1992),  in  which  images  are  represented  as  linear  combinations 
f  ^  symmetric  basis  functions.  We  have  some  preliminary  experience  with  the  behavior  of 

the  EM-ML  algonthm  using  smooth  (continuous  value  and  derivative  at  the  basis  function  “border”) 
sphencally  symmetric  functions  {blobs,  Matej  and  Lewitt  1992).  These  experiments  have  shown  that 
incorporation  of  blobs  into  the  reconstruction  model  can  be  advantageous,  especiaUy  for  the  case  of 
projection  data  with  limited  resolution  (smoothed  by  physical  processes  during  measurement).  Such 
smoothed  projection  data  have  better  consistency  with  the  reconstruction  model  incorporating  blobs. 
Moreover,  this  type  of  image  representation  can  be  advantageous  also  from  the  point  of  view  of  3-D 
visudization  algorithms  (Udupa  and  Herman  1991).  We  implemented  both  voxel-based  and  blob-based 
versions  of  EM-ML  and  of  ART  for  the  fully  3-D  case. 
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Specifically,  the  blobs  used  in  the  reconstructions  in  this  paper  are  based  on  the  generalized 
Kaiser-Bessel  window  functions  having  the  form  (Lewitt  1990) 


0  <  r  <  a  ? 
otherwise , 


where  r  is  the  radial  distance  from  the  blob’s  center,  Im  denotes  the  modified  Bessel  function  of  the 
order  m,  a  is  the  radial  dimension  of  the  blob  and  a  is  the  “taper”  parameter  controlling  the  blob 
shape.  The  following  blob  parameters  (Matej  and  Lewitt  1992)  were  used  in  our  experiments;  m 
=  2  (continuous  derivative  at  the  blob  borders),  radius  a  =  4.0  mm  (which  is  twice  the  length  of  a 
voxel  edge),  and  taper  parameter  a  -  10.4  (chosen  so  that  the  radial  frequency  spectrum  has  value 
zero  at  sampling  frequency  and  is  also  bounded  by  a  function  that  decays  particularly  fast  for  higher 
frequencies).  At  the  time  of  originally  designing  the  experiments,  we  decided  (without  considering 
alternatives)  to  place  the  blobs  in  the  three-dimensional  space  so  that  the  blob-centers  are  at  the  same 
locations  as  the  voxel-centers. 

The  EM-ML  method  starts  with  an  image  vector  x®,  all  of  whose  components  have  the  same 
positive  value,  and  then  proceeds  with  iterative  steps  given  by  the  formula: 


E 


Vi 


Ei  n  -  •  <  a*  X*  > 


where  x|  is  the  j\h  element  of  the  ^h  iterative  estimate  of  the  image  vector,  y,  is  the  data  item 
associate  with  the  ith  projection  bin,  aij  is  the  contribution  of  the  yth  basis  function  to  this  data 
item,  and  <  a‘,x*  >  is  the  inner  product  (forward  projection)  of  the  ixh  row  of  the  matrix  A  and 
the  ^  iterative  estimate  of  the  image  vector.  In  our  implementation  of  the  EM-ML  method,  the 
system  of  equations  was  set  up  so  that  y,-  is  the  actual  number  of  detected  counts  associated  with  the 
ith  projection  bin.  This  means  that  in  determining  entries  of  A  we  have  to  take  into  consideration 
the  attenuation  that  takes  place  in  the  object  whose  activity  distribution  we  wish  to  reconstruct.  This 
is  the  appropriate  model,  since  the  EM-ML  method  as  described  above  converges  to  the  maximum 
likelihood  estimator  of  the  image  vector  only  if  the  y,-  are  samples  from  Poisson  distributions  whose 
means  are  <  a*,x  >.  (By  the  same  token,  the  system  should  contain  equations  for  only  those  bins 
for  which  there  is  a  collected  measurement.)  In  practice  we  have  to  stop  the  algorithm  at  some  point; 
the  parameter  that  we  leave  free  in  our  initial  description  of  the  EM-ML  method  is  the  number  of 
iterations  at  which  the  process  is  terminated. 

The  ART  method  starts  with  an  image  vector  x®,  all  of  whose  components  have  the  same  positive 
value  (which  is  the  mean  value  over  the  reconstruction  region  as  estimated  from  the  projection  data), 
and  then  proceeds  with  iterative  steps  given  by  the  formula: 


Vi' 


<  a\  > 


1  a  ■ 


where  A  is  the  so-called  relaxation  parameter,  which  is  the  parameter  we  leave  free  at  this  initial 
stage  of  describing  this  algorithm,  and  everything  else  is  interpreted  as  it  was  done  for  EM-ML. 
(In  ART,  just  as  in  EBP,  y,  is  the  estimated  integrated  activity  associated  with  the  ith  bin;  this  is 
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obtained  from  the  detected  counts  by  a  correction  for  attenuation,  Furuie  et  al  1994.)  Dissimilarly  to 
the  EM-ML  method  (in  which  the  estimated  image  vector  is  updated  simultaneously  based  on  all  the 
measurements),  each  update  in  ART  of  the  estimated  image  vector  is  based  on  a  single  projection  bin, 
and  the  method  cycles  through  all  the  bins  one  after  another.  This  leads  to  a  faster  initial  approach 
to  ^  acceptable  solution  and  consequently  to  a  lower  number  of  necessary  total  cycles  through  the 
projection  data.  (In  the  experiments  described  in  (Matej  et  al.,  1994),  only  one  cycle  through  the 
projection  data  was  used,  except  that  it  started  with  and  repeated  at  the  end  corrections  for  the  bins 
associated  with  untilted  planes.  For  the  particular  mode  of  data  collection  that  was  simulated,  this 
means  that  the  reconstruction  process  cycles  1.07  times  through  all  of  the  projection  data.)  During 
cycling  tough  the  data  we  use  a  special  ordering  of  the  bins  (Herman  and  Meyer  1993),  which  results 
in  a  rapid  early  progress  towards  an  acceptable  solution.  Corrections  are  made  for  only  those  bins 
for  which  there  is  a  collected  measurement.  We  fixed  the  way  we  cycle  through  the  data  for  fairness 
of  comparison  between  the  various  basic  approaches:  we  wished  to  have  only  one  free  parameter  for 
which  we  optimize  during  training.  Our  experience  with  ART  indicated  to  us  that  more  is  likely  to 
be  gained  in  fixing  the  number  of  iterations  and  optimizing  for  the  relaxation  parameter  than  doing 
the  opposite.  We  expected  that  a  very  small  number  of  cycles  tough  the  data  will  give  acceptable 
results  (Herman  and  Meyer  1993). 
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Chapter  4.  IMPLEMENTATION  DETAILS 


In  this  chapter  we  describe  in  details  the  phantom  model,  reconstruction  volume  within  which  we 
do  our  reconstruction,  the  digitization  of  this  volume  into  a  voxel  space,  the  arrangement  of  lines  for 
which  projection  data  are  to  be  generated,  and  the  method  of  raysum  generation,  including  noise. 


4.1.  Reconstruction  volume 


1.  It  is  assumed  that  the  reconstruction  volume  (fig.  5)  is  a  cylinder  with: 

Dxy=256  mm  (diameter), 

Hz=128  mm  (height). 

2.  The  origin  of  a  Cartesian  coordinate  system  (xyz,  right  hand  convention)  is  located  at  the  center  of 
the  cylinder,  where  the  z-axis  coincides  with  the  cylinder  axis  (figure  5).  In  order  to  be  compatible 
with  medical  notation  for  images  (transversal,  sagittal  and  coronal  views ),  the  direction  of  z-axis  is 
from  top  to  bottom.  More  specifically,  the  patient  is  laying  on  xz-plane,  with  feet  in  +z  direction, 
looking  at  +y-axis  direction,  and  with  left  ear  in  +x  direction. 


Nt:  Num.  of  tilts 

Na:  Num.  of  angles  of  view  per  tilt 
Nr:  Num.  of  rows  per  angle  of  view 
Ns:  Num.  of  samples  per  row 


Figure  5.  Diagram  of  reconstruction  volume  (cylinder),  voxel  space  (cuboid)  and  panel  of  sensors 
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4.2.  Voxel  space 

1.  Consider  a  cuboid  that  tighdy  surrounds  the  cylinder  which  is  the  reconstruction  volume. 

2.  This  cuboid  is  divided  into  small  cubic  voxels  of  side  VoxelSide=2mm. 

3.  Intensity  range  of  values  to  be  assigned  to  voxels:  16bits  (short  integer,  2B). 

4.  Thus  we  have: 

Nx=128  (columns,  x-direction), 

Ny=128  (rows,  y-direction), 

Nz=64  (slices,  z-direction), 

and  NxxNyxNz=lM  (requiring  a  total  of  2MB  to  store). 

5.  The  voxel  at  position  (ix,iy,iz)  is,  using  C  notation: 

V[iz][iy][ix], 

where  iz=0..(Nz-l),  iy=0..(Ny-l),  and  ix=0..(Nx-l). 

6.  Correspondence  with  reconstruction  volume  (see  figure  5,  indices  ix,iy,iz): 

The  first  slice  corresponds  to  the  top  of  ’reconstruction  volume,’  i.e.,  to  the  most  negative  z  region. 
For  each  slice,  the  rows  go  in  opposite  direction  in  relation  to  y-axis  (the  first  row  starts  at  y=Dxy/2 
mm),  and  the  columns  increase  as  x  increases  (the  first  column  starts  at  x=-Dxy/2  mm).  Thus, 
center  of  voxel  V[0][0][0]  is  at  (x,y,z)=(x(0),  y(0),  z(0))  mm,  where 
x(0)=-(Dxy/2-VoxelSide/2)  mm, 
y(0)=Dxy/2-VoxelSide/2  mm, 
z(0)=-Hz/2+VoxelSide/2  mm, 
and  center  of  voxel  V[iz][iy][ix]  is  at: 
x(ix)=(ix-Nx/2+0.5)xVoxelSide  mm;  ix=0..(Nx-l), 
y(iy)=(Ny/2-0.5-iy)xVoxelSide  mm;  iy=0..(Ny-l), 
z(iz)=(iz-Nz/2+0.5)xVoxelSide  mm;  iz=0..(Nz-l). 

7.  For  display  purpose,  V[k][0][0]  is  located  at  upper  left  comer  (in  axial  view). 

4.3.  Phantoms 

1.  Each  phantom  is  contained  within  a  large  ellipsoid  with  its  z-axis  coinciding  with  the  z-axis  of 
xyz  coordinate  system  and  its  center  at  (0,0,Hz/2)  in  mm.  Between  the  planes  z=Zmin=-Hz/2 
mm  and  z=Zmax=Hz/2  mm,  there  are  nObjects=69  spheres,  distributed  in  three  layers  (A,B,C) 
(see  section  3.1).  These  spheres  contain  the  objects  (features)  that  are  ellipsoids.  The  activity  is 
restricted  to  the  reconstruction  region(fig.l). 

2.  The  radii  of  the  large  ellipsoid  are: 

Rx=128  mm, 

Ry=96  mm, 

Rz=147.80  mm. 

The  associated  intensity  is  BKG  (background).  This  value  will  be  determined  so  that  the  total 
detected  activity  (attenuated  data)  from  the  large  ellipsoid  (background  only)  will  be  80  million 
counts.  Note:  whole  phantom  is  considered  for  attenuation,  i.e,  if  a  ray  (tube)  also  intercepts  part 
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of  phantom  without  activity,  line  integral  of  activity  will  consider  only  the  part  of  phantom  with 
activity,  while  attenuation  is  computed  on  the  basis  of  the  whole  intersection  segment. 

3.  Each  feature  is  inside  a  sphere,  whose  position  and  radius  are  described  by  (xcenter,  ycenter,  zcen- 
ter,  sphereRadius)  and  whose  background  intensity  is  given  by  bkg  (local  additional  background, 
usually  equal  to  0).  See  appendix  D  for  list  of  sizes  and  locations  of  spheres. 

4.  Each  feature  (ellipsoid),  has  the  following  characteristics: 

a.  same  center  as  the  enclosing  sphere; 

b.  radii  (rx,ry,rz)  are  between  50%  to  100%  of  sphereRadius,  each  independently  randomly 
selected  from  a  uniform  distribution  over  this  range; 

c.  orientation:  the  ellipsoid  above  described  is  rotated  by  (in  this  order): 

ax  (rotation  around  x-axis), 
ay  (rotation  around  y-axis), 
az  (rotation  around  z-axis), 

where  each  of  these  three  angles  is  independently  randomly  chosen  from  a  uniform  distribution 
over  the  range  [0,7r). 

d.  activity  is  randomly  assigned  to  be  one  of  three  equiprobable  types:  NORMAL  (the  activity 
in  the  feature  is  the  same  as  in  its  background),  HOT  (the  feature  is  a  “hot  spot”),  or  COLD 
(the  feature  is  a  “cold  spot”).  The  activity  assigned  to  hot  spots  (respectively,  to  cold  spots) 
is  higher  (respectively,  lower)  than  the  activity  in  their  background  by  the  contrast  listed  in 
Table  I. 

e.  each  feature  contains  at  least  5  voxels  inside  and  5  voxels  in  feature  background  region 
(region  between  feature  and  enclosing  sphere) 

5.  In  order  to  accmately  estimate  the  average  intensity  in  a  voxel  in  the  phantom  image,  each  voxel 
is  subdivided  into  64  sub-voxels  (fig.6)  and  the  intensity  of  each  one  of  these  is  calculated  as 
the  sum  of  the  intensities  of  all  features  and  spheres  which  contain  subvoxel  center.  The  voxel 
intensity  is  the  average  of  the  intensities  of  its  sub-voxels. 

Voxel  and  sub-voxels  ((nSubVoxelsPerSlde=4) 


Figure  6.  Illustration  of  voxel  and  subvoxels 

6.  Algorithm  for  PhantomSD  creation  (.img) 

a.  Generate  (randomly)  feature  parameters  (size,  orientation,  activity) 
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b.  Check  for  at  least  5  voxels  inside  each  feature  and  background  region:  if  not,  generate 
(randomly)  new  feature.  We  are  considering  “inside”  if  center  and  at  least  five  vertices  of 
voxel  are  inside  specific  region. 

c.  For  each  feature 

•  Get  featiue  characteristics  {  rx,ry,rz,  ax,ay,az,  value  } 

•  Update  values  of  sub-voxels  inside  the  sphere  by  feat.bkg  (not  needed  if  bkg=0) 

•  Update  all  sub-voxels  inside  this  feature  by  featvalue 

d.  Calculate  each  voxel’s  value  as  average  of  the  values  of  its  sub-voxels 

e.  Store  data  (.pde  and  Petview  format  .img) 

4.4.  Projection  data 


z 


origin  O  and  with  normal  n 

Figure  7.  Coordinate  system  (u,v)  of  panel  in  relation  to  xyz-coordinate 
system  of  reconstruction  volume.  Vector  n  is  normal  to  the  panel. 

1.  It  is  assumed  that  data  are  rebinned  into  parallel  raysums,  and  activity  are  restricted  to  the  region 
between  z=[-Hz/2,H2/2]  mm.  We  are  considering  all  detectable  rays  for  the  given  scanner  geometry 
(cylinder): 

SCNdiameter=420mm, 

S  CNheight=25 6mm. 

We  are  assuming  that  the  center  of  phantom  and  reconstruction  volume  coincide  with  scanner 
center. 
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2.  There  is  a  rectangular  matrix  (of  size  nSCNrowsxnSCNcols),  called  a  panel,  of  sensors 
(SensorSidexSensorSide  each),  which  is  rotated  to  provide  views  and  tilted  to  provide  tilts  (figure 
5),  with: 

SensorSide=2  mm, 
nSCNrows=90  (number  of  rows), 
nSCNcols=128  (number  of  columns), 

nSCNviews=96  (number  of  views,  equally  spaced  between  0  and  x), 
nSCNtilts=15  (number  of  tilts), 

’nitStep=3.75  (degrees). 

For  each  tilt  and  for  each  view,  we  calculate  for  each  sensor  a  rayswn  and  the  collection  of  such 
raysums  is  referred  to  as  the  projection  data. 

nSCNrows  was  calculated  such  that  the  panel  can  detect  all  rays  detectable  by  the  scanner  for 
all  considered  tilts: 


a.  For  a  given  tilt ’t’,  for  all  views  of  this  tilt,  the  maximum  region  of  detectable  rays  is  limited 
by  the  projection  of  two  circles:  top  (bottom)  of  scanner  cylinder  (limiting  curve)  and  of 
reconstruction  cylinder  (expanding  curve)  (see  figure  8). 

b.  The  limiting  curve  for  all  rays  from  scanner  cylinder  that  can  be  detected  is  given  by: 


^cos{t)  —  sin{t).'\J  Rl  — 


c.  Expanding  curve  for  all  possible  rays  for  given  tilt  from  reconstruction  region  (cylinder)  is 
given  by: 


K 

2 


cos{t)  +  sin 


—  u 


2 


where: 

hs:  height  of  scanner  cylinder 
Rs:  diameter  of  scanner  cylinder 
hr:  height  of  reconstruction  cylinder 
Rr:  diameter  of  reconstruction  cylinder 
t:  tilting  angle  (t  >=0) 
u;  panel  abscissa 
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:  Limiting  ellipse  (Scanner) 

I  ^ 

T' 

^  ■  ^^xpanding  ellipse  (reconstr.) 


^3^  Panel  mask  for 
■  a  view  and  tilt 


\  I 

y  ' 


Figure  8.  Illustration  of  panel  mask  for  bins  with  detectable  counts  for  a  given  tilt. 


3.  It  is  assumed  that  the  center  line  of  the  panel  is  in  plane  z=0. 

4.  Rotation  of  the  panel  is  around  the  z-axis  (nSCNviews),  using  the  right-hand  rule  direction.  The 
first  panel  is  at  angle  0.  (In  case  of  no  tilting,  this  corresponds  to  the  yz  plane.) 

5.  Ulting  is  around  the  middle  line  of  the  panel  (figure  5),  with  nSCNtilts  tilt  positions. 

It  is  assumed  nSCNtilts  is  odd.  TiltStep=3.75  degrees.  So  the  tilting  varies  from 
-int(nSCNtilts/2)x'IiltStep  to  int(nSCNtilts/2)xTiltStep  and  index=0..(nSCNtilts-l).  Thus  in¬ 
dex  of  ’no  tilting’  is  int(nSCNtilts/2). 

6.  Calculation  of  raysums: 

First,  raysums  at  finer  resolution  (’s’  ,  1mm  apart)  are  calculated  as  accumulated  line  integrals 
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through  features  for  each  SensorSide^  position.  The  effective  ray  sum  for  a  sensor  ’p’(2  mm  apart) 
will  be  the  local  average  and  sub-sampling  of  filtered  raysums  ’s’  (fig.  9),i.e: 


r“ 

I — 

I  o 


I  o 


0*1--- 

J 


s,  r 

-  P 


Figure  9.  Illustration  of  panel  sampling:  raysum(p)  and  panel  sub-sampling:  sub-raysum  (s,r) 
a.  Calculate  ’s’ 

Given  a  ray  in  3D  space,  calculate  the  non-attenuated  raysum  through  the  given  phantom, 
which  is  a  set  of  ellipsoids. 


b.  Filter  ’s’  to  get  ’r’ 


n  = 


S  j 
j 


j 

where  j:  over  local  VALID  neighbourhood 


Filter  (w)  is  radially  symmetric  with  Gaussian  shape  of  size  13x13  points  and  with  FWHM=4 
mm: 


_  FWHM 
(2\/2h^^ 
u,v  =  {—6. .6mm} 


Projection  smooth  is  done  to  take  into  account  inaccuracies  in  rebinning  and  in  locating 
position  of  source  (positron). 
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c.  Local  average  and  sub-sampling  of  ’r’  yielding  ’p’  (fig.  9) 

7.  Noise  (see  section  3.2) 

8.  Store  each  panel  data  in  Petview  format  (.sen)  by  rows  p[iR],  iR=0..nSCNrows-l  (fig.  7) 
Algorithm: 

For  each  panel 

a.  Consider  sub-raysums  (SensorSide/2  by  SensorSide/2) 

b.  Reset  all  sub-raysums  to  zero 

c.  For  each  feature=0..(nObjects-l)  (including  spheres  and  large  ellipsoid),Do 

{  Find  all  rays  that  pass  through  the  feature; 

For  each  such  ray,  calculate  the  line  integral  (length  of  the  segment  defined  by  the  intersection 
of  the  feature  and  the  line,  multiplied  by  the  intensity)  through  the  feature; 

Accumulate. 

} 

d.  Filter  (Gaussian) 

e.  Calculate  raysums  as  average  of  sub-raysums  that  form  each  raysum 

f.  Noise  generation 

g.  Store  data  (Petview  format) 
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Chapter  5.  FILES  AND  FORMATS 


5.1.  Files 

1.  filename.pde:  phantom  and  scanner  geometry  description  file  (described  below) 

2.  filename.img:  3D-phantom  data  (Petview  format) 

3.  filename.scn:  3D-phantom  projection  data  (Petview  format) 

4.  phaLengths.AttenF.scn:  attenuation  factor  for  each  projection  bin  (Petview  format).  In  our  case, 
it  contains,  for  each  ray,  where  L  is  the  intersection  length  (in  mm)  of  the  ray  and  the 
phantom,  considered  as  a  non-truncated  ellipsoid  and  /x=0.(X)95  mm~^ 

5.2.  Archiving  format  for  .pde 

A  phantom  and  scanner  description  file  (.pde,  for  generation  of  masks  and  for  evaluation  purposes) 
consists  of  a  header  in  the  format  of  the  structure  PHA_HDR,  followed  by  nObjects=69  descriptions  of 
features  in  the  format  of  the  structure  FEATURE  and  structure  USER_PARAM.  Precise  specifications 
of  these  structiues  are  given  in  Appendix  C. 

5.3.  Archiving  format  for  .img  and  .sen 

The  phantom  image  data  (.img)  and  projection  data  (.sen)  are  stored  using  Petview  format  (see 
Appendix  E). 

In  order  to  make  comparison  simpler,  we  are  normalizing  background  values  to  1  “arbitrary 
property  unit”  per  mm  (relative  counts/mm  of  detector  “tube”)  and  using  a  CT  model  for  projection 
data.  Therefore  the  expected  values  for  voxels  (qj)  in  background  region  is  1  and  if  a  ray  intercepts 
256  mm  of  background  region  with  activity,  the  expected  projection  value  (pi)  is  256.  However,  we 
can  convert  back  to  realistic  units  and  to  PET  model  using  simple  factors  as  shown  below. 

The  parameter  ’TotalCounts’  provided  by  the  user  (defaults  is  80  millions)  refers  to  the  total  counts 
detected  (attenuated  projection  data)  by  a  scanner  due  to  only  background  in  phantom  .  Since  we  know 
the  geometiy  of  the  phantom  and  scanner,  and  linear  attenuation  coefficient  (/x=0.0095  mm~^),  we 
can  obtain  the  relationship  between  1  “arbitrary  property  unit”  and  counts/voxel  and  counts/bin.  This 
factor  is  stored  in  file  .pde  as  Param.BKGnew  (0.292  for  default  parameters)  and  as  pha_hdr.BKG. 


Param.BKGnew  — 


TotalCounts 


-ti.L- 


where: 

L,-:  intersection  (in  mm)  of  ray  ’i’  with  entire  phantom,  (attenuation  part) 

Li  intersection  (in  mm)  of  ray  ’i’  with  phantom  limited  by  reconstruction  region  (fig.  1)  (active 
part  of  relative  phantom). 

,• :  sum  over  all  possible  (not  missing)  rays. 
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If  we  call  : 

rj:  value  of  voxel  ’j’  [relative  counts/mm  of  “tube”] 

Pi:  value  of  projection  ’i’  [relative  counts  in  bin] 

L-j  :  length  of  intersection  of  i-th  ray  with  j-th  voxel  [mm] 

CT  reconstructin  model  tries  to  solve  : 

.L .  ..r  . 

J  J 

We  will  be  referring  to  projection  value  (pi,  yi,  gi)  as  the  value  of  projection  without  attenuation, 
or  if  one  prefers,  projection  value  after  attenuation  correction. 

If  one  wants  to  use  data  in  CT  model,  but  in  realistic  units  per  mm  of  “tube”,  it  is  enough  to 
consider: 

Xj  =  Param.BKGnew  x  rj  (for  voxel)  [density/mm  of  tube] 
yi  =  Param.BKGnew  x  pi  (for  projections)  [density  in  bin] 


Pi 

[relative  counts 
in  bin] 


y 

I 

[density  in  bin] 


For  PET  model  of  projection  data,  we  need  to  obtain  the  relationship  between  probability  for  a  bin 
’i’  to  detect  a  photon  emitted  at  location  ’j’  (tij)  and  length  of  intersection  of  ray  ’i’  and  voxel  ’j’  (Ly). 
Let  : 

fj:  number  of  photons  emitted  at  location  ’j’  [counts  in  voxel] 
gi:  number  of  photons  detected  at  bin  ’i’  [counts  in  bin] 

A:  area  of  one  detector  [mm^] 

V:  volume  of  one  voxel  [mm^] 


We  have: 


9  ■ 

t 


E 

^  J 


j 


E.t . .  =1  for  all '/  with  no  missing  data 

i  tj 
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Since  ^  ~  ^.nSCNtilts.nSCNviews  =  F  [mm/voxet\  (F=2880  for  defaults  parameters) 

we  can  approximate  t-  L-jF 

Then  Comparing  with  equation  =  Y  realistic  CT 

model  data: 

gi  =  yi 

fj  =  F  X  Xj 

Therefore,  realistic  values  of  PET  model  data  can  be  approximated  from  the  relative  CT  model 
data  by: 

fj  =  Param.BKGnew  x  F  x  rj  (for  voxel)  [counts  in  voxel] 
gi  =  Param.BKGnew  x  pj  (for  projections)  [counts  in  bin] 


9i 


[counts  in 
bin] 


Approximate  PET  model 


CT  model 
Reconstruction 


! 


f . 

I 


[counts  in 
voxel] 


Some  possible  scenarios  for  using  our  projection  data  and  reconstruction  programs: 

1.  Using  our  projection  data  pi  (relative  CT  data  model)  and  provided  reconstruction  programs  (for 
CT  model) 

a.  We  will  get  reconstructed  values  rj  in  relative  values,  i.e.,  the  expected  value  in  background 
is  1. 

b.  If  we  want  results  in  “density  per  mm  of  tube”  we  should  multiply  input  or  output  data  by 
Param.BKGnew,  we  get: 

Xj  =  Param.BKGnew  x  rj 

c.  If  we  want  results  in  counts  per  voxel,  we  should  multiply  input  or  output  data  by 
Param.BKGnew  x  F,  we  get: 

fj  =  Param.BKGnew  x  F  x  rj 

(that  is  approximately  what  we  would  get  if  we  have  used  a  PET  reconstruction  model  with 
projection  data  :  gi  =  Param.BKGnew  x  pi) 
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2.  Using  our  projection  data  (pi  )  associated  to  a  phantom  (rj  )  with  other  reconstruction  programs: 

a.  If  it  uses  a  CT  model  for  reconstruction,  it  is  enough  converting  our  projection  data  (Petview 
format)  to  program’s  format  (the  result  will  be  in  relative  units,  rj). 

b.  If  it  uses  a  PET  model  for  reconstruction,  the  reconstructed  image  will  be:  F  x  rj  [relative 
counts/voxel].  In  order  to  obtain  proper  number  of  counts/voxel,  we  should  multiply  results 
by  Param.BKGnew. 

3.  Using  provided  reconstruction  programs  (CT  model)  with  other  projection  data: 

a.  CT  model  projection  data: 

•  set  Param.BKGnew=l 

•  can  be  used  directly  and  the  background  value  of  reconstructed  image  will  be  actual 
background  value  [density/mm  (of  tube)]. 

b.  PET  model  for  projection  data  (gi)  associated  to  an  image  (fj)  [counts  per  voxel]: 

•  set  Param.BKGnew=l 

•  the  reconstructed  image  will  be  :  ^ 

Notes: 

•  In  the  noise  generation  discussed  in  previous  sections,  we  used  values  of  ’gi’. 

•  Our  data  are  stored  in  CT  model  — voxel  values  (rj)  and  projection  data  (pi).  Before  storing  they  are 
scaled  (multiplied)  by  ’Param.Image_scale’  and  ’Param.MultFactorForStore’,  respectively.  Default 
values  of  scaling  factors  are  1(XX).  This  procedure  was  necessary  to  overcome  the  equalization 
errors  of  transforming  real  number  into  integer,  since  Petview  format  (details  in  pet_mhdr  and 
pet_sbhdr  structures  in  file  pet_hdr.h.  Appendix  E)  for  data  is  integer  (2  bytes).  In  Petview 
write  functions,  the  numbers  in  each  slice  are  automatically  scaled  (by  sh.imgscl  and  sh.scnscl, 
respectively)  to  not  overflow  integer  type. 

More  specifically,  archiving  is  done: 

1.  Phantom  (.img) 

a.  Storage  is  slice-oriented,  usually  from  top  to  bottom  of  structure. 

b.  For  each  slice  (z  fixed),  the  data  are  stored  row  by  row,  i.e., 

for(iy=();  iy<Ny;  iy-l-+) 

for(ix=0;  ix<Nx;  ix-H-)  store  r[iz][iy][ix]  x  Param.Image_scale 
(Note:  function  wrimgO  automatically  calculates  sh.imgscl  and  scales  for  each  slice  separately. 
So  the  stored  value  in  file  is  r[iz][iy][ix]  x  Param.Image_scale/sh.imgscl,  type  int,  where  sh 
is  a  struct  of  type  pet_sbhdr).  However,  Petview  function  rdimgO  does  not  automatically 
rescale  the  read  data. 

Size=Nz*Ny*Nx=lM  (2MB) 

2.  Projection  data  (.sen): 

a  frame  parameter  is  a  set  of  panels  for  a  particular  tilt  (there  are  nSCNtilts  frames), 
each  frame  has  nSCNviews  x  nSCNrowsxnSCNcols  raysums; 
a  slice  parameter  is  a  set  of  raysums  for  a  particular  view  (panel),  for  a  particular  tilt  (there  are 
nSCNviews  “slices”  per  frame). 
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each  slice  has  nSCNrows  x  nSCNcols  raysums; 

C  notation: 
p[iT][iV][iR][iC] 

iT=0..(nSCNtilts-l); 
iV=0..(nSCNviews-l); 
iR=0.  .(nS  CNrows- 1 ); 
iC=0..(nSCNcols-l). 

Algorithm: 

For  each  frame,  store  slice  by  slice, 

For  each  slice  (panel),  store  row  by  row,i.e.: 
store  p[iT][iV][iR][iC]  x  Param.MultFactorForStore 
where  row  start  at  most  negative  v  (iR=0)  (figure  7) 
and  column  start  at  most  negative  u  (iC=0). 

(Note:  function  wrimgO  automatically  calculates  sh.scnscl  and  scales  for  each  slice  separately.  So 
the  stored  value  in  file  is  p[iT][iV][iR][iC]  x  Param.MultFactorForStore/sh.scnscl,  type  int,  where 
sh  is  a  struct  of  type  pet_sbhdr).  However,  Petview  function  rdimgO  does  not  automatically 
rescale  the  read  data. 

Size=nSCNtilts  x  nSCNrows x  nSCN views  x  nSCNcols=  1 6.6M  (33.2MB) 
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Chapter  6.  FLOW  DIAGRAM  for  genPhantomSd.c 


.Generate  noise  (raysums) 
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Appendix  A  Justification  for  contrast  function 

A  common  index  for  measuring  detectability,  when  frequency  distribution  for  a  given  decision 
variable  is  nearly  Gaussian,  is  detection  signal-to~noise  ratio  given  by  (Fukunaga,  1972;  Hanson,  1990): 


where  ctq  ,  are  average  and  variance  of  signals  xi  (feature)  andxo  (background)  respectively. 

For  reconstructed  images,  it  can  be  assumed  that  the  noise  variance  of  voxel  intensities  inside 
features  and  outside  are  equal  and  independent  of  the  size  of  features.  Then, 


=  fi _ 

cr 

where  cr^  =  <^i  =  (coiistant  over  whole  phantom). 
For  a  given  feature  f  with  contrast  c/,  we  have: 


Then, 


E{xi)  E(xg)  _  •£'(2:1)  E(xq)  _  Cy 


Var  d  .  = 


Var(xj)  +  yar(a:Q) 


Wq  + 

ngMj 


where  we  assumed  that  the  noise  variance  is  constant  and  nQ,nj  are  number  of 

voxels  considered  for  averaging. 

In  the  analysis  of  detectability,  for  a  given  reconstruction,  we  are  using  phantom  features  with 


h  ^ 

different  sizes.  Therefore,  it  is  important  to  have  approximately  the  same  detection  S/N  ratio 

for  all  the  featiue  sizes.  This  can  be  accomplished  by  setting  the  contrast  of  features  proportion^  to 
their  sizes: 


P+2X  (23) 

/  V  "o"! 

For  our  case,  the  feature  is  inside  a  sphere  that  provides  background  region  for  the  featiue.  Thus, 
the  number  of  voxels  is  proportional  to  R  and  consequently  Cr  (x  R  2 .  However  due  to  other 

_ 2  , 

effects  (such  as  border  effect,  randomness  of  sizes)  experimentally  we  found  that  c  ^  cx  R  yields 
comparable  detectability  for  different  sizes.  In  our  experiment  we  considered  3  ranges  of  enclosing 
spheres  sizes  ([4.0  mm,  6.4  mm],  [8.0  mm,10.0mm],  [16.0  mm,  20.0mm])  and  we  used  the  function 

c  - 

""f  -  W 


36 


Appendix  B  Justification  for  detectability  index  formulae 


In  Appendix  A  we  have  derived  a  theoretical  fOTmula  for  the  feature  contrast,  enabling  comparable 
detectability  for  features  of  different  sizes.  Since  we  are  going  to  use  ROC  analysis  with  features  with 
different  sizes,  it  is  important  to  have  also  and  Var(^dj'j  approximately  independent  of  the 

size.  It  can  be  achieved  by  scaling  detectability  index  df  by  contrast  function.  . 

A  new  detectability  index  can  be  defined  by; 


h 


f 


Considering  contrast  function  (Appendix  A)  as: 


c^=k. 


f»o 

n^.nj 


Then  E{h^^  and  Var(ji  will  be  independent  of  the  feature  size: 


Var 


We  have  (k=l): 


h 


f 


Var 


o-ni 


+  n 


1 


=  1 

(^1  ~  ^o) 
a 


One  estimate  of  a  is  (Meyer,  1975): 

^0  +  *^1 

where  =  ^.  E  (r^oi  “  '  “d  o?  =  ^-  E  (r^l,  -  ^7)^ 

Hence  we  have  equation  (4): 


h 


f 


(^1  ^o) 


\f^i 


.V 


2 

1 
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Appendix  C  Definitions  (flle:genPhantom3d.inc) 


File:  genPhantom3d.inc 

/*— - 

Definitions  for  reconstruction  volume  and  voxel  space 

- - - — - */ 

#define  _Dxy  256  /*  mm,  diameter  of  reconstr.cylinder  */ 

#define  _Hz  128  /*  mm,  height  of  reconstr.  cylinder  */ 

#define  _VoxelSide  2  /♦  mm  ♦/ 

#define  _nSubVoxelsPerSide  4  /♦  64  sub-voxek  pa  Voxel  */ 

#define  _nObjccts  69 

#define  _Nx  LDxy+_VoxelSide-l)/_VoxelSide 
tdefine  _Ny  _Nx 

#define  _Nz  CHz+_VoxelSide-l)/_VoxelSide 

/* - - - 

Definitions  for  truncated  ellipsoid  and  features 

- - - - - */ 

#define  _Zmin  -  _Ez/2 

#define  _Zmax  _H2/2 

#define  _Rx  128.0  /*  mm,  radius  in  x-axis  */ 

#define  _Ry  96.0 
#define  _Rz  147.80 

#define  _BKG  1.00  /*  Initial  relative  value  .Later  it  is  set  equal  to  Param.BKGnew*/ 
#define  _hA  60.0  /*  mm,  height  of  region  A  */ 

#define  _hB  40.0 
#define  _hC  28.0 

#define  _RA  20.0  /*  mm,  radius  of  large  spheres  in  A  */ 

#define  _RB  16.0  /*  mm,  radius  of  large  spheres  in  B  */ 

#define  _RC  10.0  /*  mm,  radius  of  large  spheres  in  C  */ 

#define  _rA  8.0  /*  mm,  radius  of  small  spheres  in  A  */ 

#define  _rB  6.4  /*  mm,  radius  of  small  spheres  in  B  */ 

#define  _rC  4.0  /*  mm,  radius  of  small  spheres  in  C  */ 

#define  _MinNum Voxels  5  /*  minimum  feature  size  */ 

#define  _MinNum Vertices  5  /*  minimum  feature  size  */ 

#define  _MaxNumLoops  300  /*  for  randomizing  parameters  */ 

Definitions  for  projection  data 

- - - - */ 

#define  _SCNdiameter  420  /*  mm  */ 

#define  _SCNheight  256  /*  mm  */ 

#define  _SensorSide  2  /*  mm  */ 

#define  _nSCNrows  90 

#define  _nSCNcols  128  /*  num.of  sensor  cols  in  the  scanner  */ 

#define  _nSCNviews  96  /*  num.of  views  around  the  phantom  [0..pi)  */ 

#define  _nSCNtilts  15  /♦  num.of  tilts  */ 

#define  _TiltStep  360.0/_nSCNviews  /•  3.75  degrees  */ 

#define  _TotalCounts  80.0e6  /*  total  num.of  emissions  detected  */ 

#define  _nSubSensorPerSide  2  /*  4  sub-projections  per  projection  */ 

#define  _ATTENcoef  0.0095  /*  attenuation  coef.  in  mm-1  */ 

#define  _coefProjPoisson  0.0  /*  coef.influence  of  raysum  onto  pseudoPoisson*/ 
#define  _FLTdiameter  13.0  /*  mm,  filter  diameter  ♦/ 

#define  _FWHM  4.0  /*  mm,  of  filter  */ 

#define  _MultFactorForStore  1000.0  /*  multiplier  for  storing  .sen*/ 
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#define  _Image_scale  1000.0  /*  multiplier  for  storing  .img  */ 
#(lefine  _MissingData  1  /*  true  missing  data  */ 

/* - 

Types  and  structure  definitions 

- */ 

typedef 

struct  PHA_HDR 

{ 

/♦ - 

Header  (256  B) 


- ^ - */ 

short  nObjwts;  /♦  #  of  objects 

I*  body:  truncated  ellipsoid.  Center  at  (0,0, -Hz/2)  */ 

float  Rx,RyJlz;  /Tadii  in  mm  */ 

float  BKG;  /♦  intensity  inside  truncated  ellips.*/ 

float  ATTENcoef;  /*  attenuation  of  main  body  (background*/ 

/*  Reconstruction  volume  */ 

float  Zmin;  /*  top  of  phantom,usually  -Hz/2  mm  */ 

float  Zmax;  /*  bottom  of  phantom,  usually  Hz/2  mm  */ 

float  Xmin,  Xmax; 

float  Ymin,  Ymax; 

short  Nx,NyJ^z;  /*  num.  of  elements  in  each  direction  */ 
char  comments[204  ]; 

)  pha_hdr  ; 

typedef  struct  FEATURE 

( 

/* - 

Data  (16*4=64B) 

- ^ - */ 

/*  **  object  **  */ 
float  sphereRadius; 
float  xcenter,ycenter,  zcenter; 
float  bkg;  /*  aditional  intensity  inside  sphere  */ 
float  rx,ryjz;  /*  radii  of  ellipsoid  */ 
float  ax,ay3z;  /*  rotation  of  feature  (radians) 
in  relation  to  x,y,z  axis*/ 
float  value;  /*  aditional  intensity  inside  feature*/ 
float  AttenCoef;  /*  atten.coef.of  feature  */ 
float  spare0,sparel,spare2;  /*  future  use  */ 

)  feature; 

typedef  struct  USER_PARAM 

{ 

char  prefix[30]; 

BOOLEAN  RANDOM; 
int  nPhantoms; 

BOOLEAN  POISSON; 

BOOLEAN  FOV; 

BOOLEAN  INTERSECTIONS; 
char  fileLengths[80]; 
float  coefProjPoisson; 
float  TotalCounts; 
char  comments  [200]; 
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/*  for  reconstruction  volume  and  voxel  space  */ 
float  Dxy, 

Hz, 

Voxels  ide; 

int  nSubVoxelsPerSide, 
nObjects, 

Nx, 

Ny, 

Nz; 

float  Zmin, 

Zmax; 

/*  for  truncated  ellipsoid  and  features 
BKGnew:  actual  value  of  background 
*1 

float  Rx, 

Ry. 

Rz, 

BKG,  BKGnew, 

ATTENcoef, 

hA, 

hB, 

hC, 

RA, 

RB, 

RC, 
rA, 
rB, 
rC; 

int  MinNum Voxels; 

/*  for  projection  data  */ 
float  SensorSide; 
int  SCNdiameter, 

SCNheight, 
nSCNrows, 
nSCNcols, 
nSCNviews, 
nSCNtilts, 
nSubSensorPerSide; 
float  TiltStep; 
float  FWHM; 
float  FLTdiameter; 
float  MultFactorForStore; 
float  Image_scale; 
int  MissingData; 

char  Comame[80];  /*  name  of  correction  file*/ 

/*  area  under  program  control(temporary  mem)  */ 

char  phantomN[4];  /*  phantom  #  0..999  */ 

float  factorSumRl; 

long  seed;  /*  for  numb.random*/ 

int  TESTING;  /*  flag  for  test  */ 

)  UserParam; 


Appendix  D  Table  of  sphere  sizes  and  locations 


Feature 

Sphere  radius 

(x,y.z)  of  sphere  center 

(mm) 

(mm) 

1 

20.0 

(-40.00,  40.00,  34.00) 

2 

20.0 

( 0.00,  40.00.  34.00) 

3 

20.0 

( 40.00,  40.00,  34.00) 

4 

20.0 

(-80.00,  0.00,  34.00) 

5 

20.0 

(-40.00.  0.00,  34.00) 

6 

20.0 

( 0.00,  0.00,  34.00) 

7 

20.0 

( 40.00,  0.00,  34.00) 

8 

20.0 

(  80.00,  0.00,  34.00) 

9 

20.0 

(-40.00,-40.00,  34,00) 

10 

20.0 

(  0.00,-40.00,  34.00) 

11 

20.0 

(  40.00,-40.00,  34.00) 

12 

8.0 

(-60.00,  20.00,  50.00) 

13 

8.0 

(  60.00,  20.00,  50.00) 

14 

8.0 

(-60.00,-20.00,  50.00) 

15 

8.0 

(  60.00.-20.00,  50.00) 

16 

8.0 

(-60.00,  20.00,  34.00) 

17 

8.0 

(  60.00,  20.00,  34.00) 

18  1 

8.0 

(-60.00,-20.00,  34.00) 

19 

8.0 

(  60.00,-20.00,  34.00) 

20 

8.0 

(-60.00.20:00,  18.00) 

21 

8.0 

(  60-00,  20.00,  18.00) 

22 

8.0 

(-60.00,-20.00,  18.00) 

23 

8.0 

(  60.00,-20.00,  18.00) 

24 

16.0 

(-32.00,  32.00,-16.00) 

25 

16.0 

(  0.00,  32.00,-16.00) 

26 

16.0 

(  32.00,  32,00,-16.00) 

27 

16.0 

(-64.00,  0.00,-16.00) 

28 

16.0 

(-32,00,  0.00,-16.00) 

29 

16.0 

(  0.00.  0.00,-16.00) 

30 

16.0 

(  32.00,  0.00,-16.00) 

31 

16.0 

(  64.00.  0.00,-16.00) 

32 

16.0 

(-32.00,-32.00,-16.00) 

33 

16.0 

(  0.00,-32.00,-16.00) 

34 

16.0 

(  32.00,-3X00.-16.00) 

35 

6.4 

(-48.00,  16.00,  -3.20) 

36 

( 48.00,  16.00,  -3.20) 
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(-48.00.46.00,  -3.20) 


( 48.00.-16.00,  -3.20) 


(-48.00,  16.00,-16.00) 


(  48.00,  16.00,-16.00) 


(-48.00,-16.00,-16.00) 


( 48.00,-16.00,-16.00) 


42 


Appendix  E  Petview  format 


The  data  files  for  image  and  projection  data  are  stored  as  .img  and  .sen  files,  using  Petview  format, 
whose  header  (petjnhdr)  and  subheader  (pet_sbhdr)  structure  are  described  below. 

In  order  to  read  and  write  files  in  Petview  format,  the  following  routines  are  needed,  which  are 
include  in  file  imagiol.l.c  ( imagio.h).  Please  see  documentation  in  these  files  for  further  information. 
The  included  reconstruction  algorithms  may  also  be  useful  as  examples  for  reading/writing  data  in 
Petview  format. 

Routines  for  reading: 

1.  fileopnQ  :  opens  file 

2.  getmhdrO  ;  reads  main  header 

3.  getsbhdrO  :  reads  sub-header 

4.  rdimgO  :  reads  data  (1  slice) 

5.  filcloseO  :  closes  file 
Routines  for  writing: 

1.  filcreO  :  creates  file 

2.  wrmhdrO  :  writes  main  header 

3.  wrimgO  :  writes  data  (1  slice) 

4.  filcloseO  :  closes  file 

The  following  routines  that  are  included  in  file  MIPGphaFunc.c  may  be  useful  : 

1.  mh_mitO  :  initializes  main  header 

2.  sh_initO  :  initializes  sub-header 

3.  mh_shIMGupdate():  updates  main  header  and  sub_header  structure  of  image  data  based  on 
PARAM  structure  content 

4.  mh_shSCNupdateO.  updates  main  header  and  sub_header  structure  of  projection  data  based  on 
PARAM  structure  content 

/* . File  :  pet_hdr.h . 

Piupose;  Contains  main  header  and  sub  header  definitions.  It  should  be  used  as  an  include  file  in  programs  which  use 
the  imagio  routines  to  readi/write  the  headers. 

By:  PetviewAJGM  group 
Date:  Revised  July  9,  1991 
*1 

/**  Cl  -  Comment  #1  -  eliminate  the  setting  of  these  parms  by  UNIX  localtimeO  call  in  wrmhdrO.  The  application  should 
set  these  values  only,  not  wrmhdrO. 

C2  -  replacing  the  ’ybias’  parameter  name  with  ’sclsep’.  Words  157-158  UGM  scanner  physical  slice  separation  is  1.0 
mm,  but  for  older  University  scanner,  the  separation  is  2.0  mm.  Currently:  if  (’ybias*  =  1.0)  data  acquired  on  UGM 
scanner  else  data  acquired  on  University  scanner 
C3  -  sentype  (scan  type)  contains  all  possible  .img  types  and  .sen  types 

**l 

struct  pet_mhdr  { 

short  filtyp;  /*  File  type.  1=.SCN,  2=.IMG,  3=.VID  4=.OTHER  */  /*Cl*/ 
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short  daycre;  /*  Day  of  file  creation.*/ 

/*C1*/  short  mocre;  /*  Month  of  file  creation.*/ 

/*C1*/  short  yrcre;  /*  Year  of  file  creation.*/ 

/*C1*/  short  hrcre;  /*  Hour  of  file  creation.*/ 

/*C1*/  short  mincre;  /*  Minute  of  file  creation.*/ 

/*C1*/  short  seccre;  /*  Second  of  file  creation.*/ 
short  pattyp:  /*  Scan  type,  0=test,  l=patient  */ 
short  shdtyp;  /*  SubheadCT  type,  0=UCLA/CTI,  1=UGM.  */ 
short  duiatn;  /*  Duration  of  scan  in  seconds.  */ 

/*C3*/  short  scntyp;  /*  0=not  a  scan,  l=blnk,  2=trmiss,  3=emiss  4=nonn  5=obl,  6=bsi  */ 
short  numray;  /*  Numbw  of  rays  per  projection.  */ 
short  niunang;  /*  Number  of  angles  in  collected  data  */ 
short  slcthk;  /*  Slice  thickness  in  mm.  */ 

short  isotop:  /♦  Isotope  (1=F18,  2=015,  3=C11  4=Ga-68,  5=N-13,  6=Rb-82,  7=Cu-62)  */ 

float  slq)e;  /*  Calibration  slope,  nciAnl/cts/pix/min.  */ 

float  intcpt;  /*  Calibration  intercept,  nci/ml.  default=0  */ 

short  injtim;  /*  Injection  time  offset  for  scan  start  */ 

short  nslice;  /*  Number  of  slices  per  frame.  */ 

short  nframe;  /*  Number  of  frames  collected.  */ 


char  bthday[7];  /*  Subject  birthdate,  (MMDDYY).*/ 

char  ssn[10]:  /*  Subject  social  security  number.*/ 

short  petnum;  /*  PET  number.*/ 

float  dmax;  /*  Diameter  of  region  being  imaged.*/ 

float  angmax;  /*  Maximum  acceptance  angle.*/ 

float  xO;  /*  X-coord.  of  center  of  image.*/ 

float  yO;  /*  Y-coord.  of  center  of  image.*/ 

float  zO;  /*  Z-coord.  of  center  of  image.*/ 

float  nevent;  /*  Number  of  raw  events  processed.*/ 

float  nsino;  /*  Number  of  events  stored  in  sinogram.*/ 

short  eglobjow;  /*  Global  energy  -  upper  threshold  */ 

short  eglob_up;  /*  Global  energy  -  upper  threshold  */ 

short  eloc_low;  /*  Local  energy  -  lower  threshold  */ 

short  eloc_up;  /*  Local  energy  -  upper  threshold  */ 

short  orientation;  /*  0  -  not  applicable,  l=Head  First,  2=Feet  First*/ 

char  scan_swrBl[6];  /*  scanner  software  release;  range  01.00-99.99  */ 

/*C2*/  float  slcsep;  /*  Physical  Slice  Separation  1  =  UGM  */ 
char  fctrfil[20];  /*  Factor  file  name.  (PMT  gains).*/  /*  was  16  */ 
char  baselin[20]:  /*  Baseline  file  name.  (DC  oflfsets).*/  /*  was  16  */ 
char  dstpkfl[20];  /*  Distortion  peak  file  name.*/ 
char  aqprotocol_name[20];  /*  acquisition  protocol  name  */ 

short  aqprotocoLtype;  /*  l=Emiss-Static,  2=Emiss-Dynamic,  3=Trans,  4=Gated  Cardiac,  5=Whole  Body  */ 

char  patient_name[30];  /*  patient  name  */ 

float  reslice_angl;  /*  Reslicing  (OBL)  angle  1  */ 

float  reslice_ang2;  /*  Reslicing  (OBL)  angle  2  */ 

float  reslice_ang3;  /*  Reslicing  (OBL)  angle  3  */ 

short  minslc;  /*  Minimum  Slice  number  */ 

short  maxslc;  /*  Maximum  Slice  number  */ 

short  minfinm;  /*  Minimum  Frame  number  */ 

short  maxfrm;  /*  Maximum  Frame  number  */ 

short  trailexists;  /*  1=  trailer  exists,  0  =  no  trailer  */ 

long  trailbeg;  /*  unsigned  32  bit  number  =  #  of  bytes  from  file  beginning  indicating  where  the  trailer  begins 
);  /*  end  pet_mhdr  */ 


struct  pet_sbhdr  { 

/*  currently  wrsbhdr  does  not  write  the  following  subheader  parameter:  gatint  */ 


*/ 
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short  datype;  /♦  Data  type  (l=byte,2=I*2,4=R*4orI*4,8=R*8)  */ 

short  xdim;  /*  X  dimension  (128  for  .IMG,  256  or  128  for  .SCN)  */ 

short  ydim;  /*  Y  dimension  (128  for  .IMG,  192  or  96  for  .SCN)  */ 

short  slcnum;  /*  Slice  number.*/ 

short  gatint;  /*  Gating  interval  (msec).  */ 

float  magfac;  /♦  Reconstruction  magnification  (nominally  1)  */ 

/*  magfac  used  in  analyze  anz2u.c  in  imagio.c  Convert  UCLA  to  Analyze  Format  ♦/ 

float  imgscl;  /*  Image  scale  factor  (nominally  1.)  */ 

int  imgmin;  /*  Minimum  value  in  image;  scaled  to  short  by  wrshdrO*/ 

int  imgmax;  /*  Maximum  value  in  image;  scaled  to  short  by  wrshdrO*/ 

float  scnscl;  /*  Sinogram  scale  factor  */ 

short  strhr;  /*  Starting  time  of  this  frame:  hour.*/ 

short  strmin;  /*  Starting  time  of  this  frame:  minute.*/ 

short  strsec;  /*  Starting  time  of  this  frame:  second.*/ 

int  scnmin;  /*  Minimum  value  in  sinogram;  scaled  to  short  by  wrshdrQ*/ 

int  scnmax;  /*  Maximum  value  in  sinogram;  scaled  to  short  by  wrshdrQ*/ 

short  endhr;  /*  Ending  time  of  this  frame:  hour.*/ 

short  endmin;  /*  Ending  time  of  this  frame:  minute.*/ 

short  endsec;  /*  Ending  time  of  this  frame:  second,*/ 

long  midtim;  /*  Time  from  start  of  scan  to  frame  midtime  of  frame.  */ 

/*  midtim  used  by  cmrglu.c  */ 

short  scnlen;  /*  Scan  duration  in  seconds.*/ 

float  imgsum;  /*  Total  number  of  events  in  image.*/ 

float  scnsum;  /*  Total  number  of  events  in  sinogram.*/ 

/*  additional  parameters  to  add  to  the  subheader  */ 

/*  Space  available  from  words  113  -  256  */ 
float  bgsngrt;  /*  Beginning  singles  rate  -  all  detectors  */ 
float  bgcoincrt;  /*  Beginning  coincidence  rate  -  all  bankpairs  */ 
float  endsngrt;  /*  Ending  singles  rate  -  all  detectors  */ 
float  endcoincrt;/*  Ending  coincidence  rate  -  all  bankpairs  */ 
float  deadtimefac;/*  Deadtime  Correction  Factor  default=LO  */ 
short  bedpos;  /*  Bed  position  for  whole  body  scanning  */ 

RECONINFO  reconinfo; 

);  /*  end  pet_Sbhdr  */  /*********************************J(e**>|t*>|c*J|c3|c3|e*j(c*%*>|ej|«>|c**>|e******>k**Hc*H«5|c5|e*Jfe*^ 
typedef  struct  reconinfo  ( 

/****  Reconstruction  Parameters  *****/ 

/*  General  */ 

char  recon_swrel[6];  /*  reconstruction  software  release;  range  01.00-99.99  */ 

char  analy_swrel[6];  /*  display/analysis(petview)  software  release;  range  01.00-99.99  */ 

char  recprotocol_name[20];  /*  reconstruction  protocol  name  */ 

char  insinofiie[20];  /*  sinogram  file  to  be  reconstructed  */ 

short  slc_add;  /*  slice  addition  l=Yes  0=No  */ 

short  slc_space;  /*  if  slice  add,  Slice  Spacing  in  slices  */ 

short  slc_thick;  /*  if  slice  add,  Slice  Thickness  in  slices  */ 

short  frame_add;  /*  frame  addition  l=Yes  0=No  */ 

short  frame_space;  /*  if  frame  add,  Frame  Spacing  in  frames  */ 

short  frame_thick;  /*  if  frame  add.  Frame  Thickness  in  frames  */ 

/*  Filtered  Baclq)rojection  */ 

short  fltr_type;  /*  0=NONE,  1=:HANNING,  2=GAUSSIAN  */ 
short  smoth;  /*  Smooth  Factor  for  Hanning  Filter  */ 

/*  Background  Subtraction  */ 

short  bgsub^type;  /♦  0=NONE,  l=UNIFORM,  2=NONUNIFORM  */ 
short  edge_exp;  /*  Number  of  bins  to  expand  edge  */ 
short  ang_avg;  /*  Number  of  angles  to  average  */ 
float  bck_coeff;  /*  Non-uniform  background  coefiScient  */ 
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short  bck_wid;  /*  Number  of  bins  for  fitting  tails  radially  */ 

/*  Attenuation  Correction  */ 

short  attncOT.type;  /*  0=NONE,  1=TRAN,  2=ELLIPSE,  3=READ  ♦/ 
float  atm_coef;  /*  Attenuation  Coefficient  in  lycm.  */ 
char  regfile[20];  /*  Attenuation  region  file  name  */ 

char  proc_transinofile[20]:  /*  Processed  Transmission  sinogram  file  of  attenuation  coefficients  */ 
float  skull_comp;  /♦  Ad  hoc  correction  for  effect  of  skull  ’"/ 

/*  Normalization  */ 

short  norm_type;  /*  0=NONE.  1=AXIAL,  2=EFFICIENCY,  3=BOTH  */ 
short  smp_norm;  /*  Sampling  Pattern  Removal  0=No  l=Yes  */ 
char  axnfile[20];  /*  Axial  Normalization  file  name  */ 
char  efffiormfile[20];  /*  Efficiency  Normalization  file  name  */ 

/*  Compensation  */ 

short  g^_comp;  /*  g^  compensation  l=Yes  0=No  */ 
short  g^_expand:  /*  NumbCT  of  bins  to  expand  gaps.  */ 
short  num_itCT;  /*  NumbCT  of  iterations.  */ 

short  g^_dir;  /*  Gap  insertion  direction.  (l=HORIZONTAL  2=VERTICAL).  */ 

short  g^_edge:  /*  Use  edge  information  (1=TRUE  0=FALSE).*/ 

short  g^_fthr;  /*  >Mdth  of  sinogram  gap  feathering.  */ 

short  mod_iter;  /*  Modified  last  iteration.  (1=TRUE  0=FALSE)*/ 

float  freq_mask_fwhm;  /*  FWHM  of  fireq  mask  smoothing  filter.  */ 

/♦  Decay  Correction:  */ 

short  decay_corr,  /*  Decay  correction  (1=TRUE  0=FALSE).  ♦/ 

/♦  Transmission  scan  processing.  */ 

char  transinofile[20];/*  Transmission  scan  sinogram.  */ 

char  blnksinofile[20];/*  Blank  scan  sinogram.  */ 

float  tran_ray_fwhm;/*  FWHM  of  Gaussian  used  for  transverse  smoothing.*/ 
float  tran_axl_fwhm;/*  FWHM  of  Gaussian  used  for  axial  smoothing.*/ 
short  surv_mask;  /*  Number  of  rays  to  mask  survival  prob  to  1.*/ 

)  RECONINFO; 
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Appendix  F  Script  sampies 


File:  buildalhbat 

#!  /bin/csh 

#  file:  buildall.bat  Sergio  Furuie  MIPG  22Sep93 

# 

#  This  script  builds  all  programs  related  to  evaluation  of  3d  PET 

#  It  should  be  run  in  subdirectory  ’bin* 

# 

#  Usage  for  example  for  spare: 

#  buildall.bat  spare 

#  (In  this  case,  it  creates  a  subdirectory  called  ’spare’) 

#  It  assumes: 

#  1)  environment  variable  ’dir3d’  set.  Example:  setenv  dir3d  /usr/EVAL3DPET 

#  2)  Sources  have  been  copied  to  subdirectories  of  ’sources’: 

#  ’bib*,  ’gen*,  ’rec’,  ’eval’,  ’util* 

#  Machine  dependent 
setenv  dirMach  $1 
mkdir  SdirMach 

#  Machine  independent. 

set  sourcesDir=$dir3d/sources 
set  genDir=$sourcesDir/gen 
set  bibDir=$sourcesDir/bib 
set  recE>ir=$sourcesDir/rec 
set  evaIDii^$sourcesDir/eval 
set  trainDir=$sourcesDir/train 
set  utilDir=$sourcesDir/util 

#  creating  library 
cd  bibDir 

make  -f  bib.mak 
Vm  *.0 

#  Phantom  generation 
cd  genUir 

make  -f  gen.mak 
Vm  *.o 

#  Reconstruction  programs 
cd  recDir 

make  -f  rec.mak 
^rm  *.o 

#  Evaluation  programs 
cd  evalDir 

make  -f  eval.mak 
Vm  *.o 

#  Training  programs 
cd  trainDir 

make  -f  train.mak 
Vm  *.o 

#  Utility  programs 
cd  utilUir 

make  -f  util.mak 
Vm  *.o 


#  File:  recon.bat  MIPG  /  SF  Jan94 

#  This  batch: 

#  Reconstruct  several  phantoms  using  4  methods 
unalias  cp 

#/etc/renice  10  -u  furuie 
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#  Reconstruct 

echo  "Reconstruction  log"  >  MIPGrecon.log 
foreach  i  (010  011  012  013  014  015  ) 
echo  "Processing  case""$i" 
echo  "Processing  ART  blob" 

artblob  -b  -i  phaMIPG"$i".scn  -f  test  -s  2  -n  1  -1  0.0055  -r  2  -m  2  -a  10.4 

cat  testhist  »  MIPGrecon.log 

echo  ‘hostname*  »  MIPGrecon.log 

complmgs  phaMIPG"$i".img  testimg  »  MIPGrecon.log 

cp  testimg  MIPG"$i"ab.img 

echo  "Processing  ART  voxel" 

artvox  -b  -i  phaMIPG"$i".scn  -f  test  -s  2  -n  1  -1  0.0049375 

cat  testhist  »  MIPGrecon.log 

echo  ‘hostname*  »  MIPGrecon.log 

complmgs  phaMIPG"$i".img  testimg  »  MIPGrecon.log 

cp  testimg  MIPG"$i"av.img 

echo  "Processing  EM  blob" 

emblob  -b  -i  phaMIPG"$i".scn  -f  test  -s  2  -n  7  -r  2  -m  2  -a  10.4  -t  phaLengths.AttenF.scn 

cat  testhist  »  MIPGrecon.log 

echo  ‘hostname*  »  MIPGrecon.log 

complmgs  phaMIPG"$i".img  testimg  »  MIPGrecon.log 

cp  testimg  MIPG"$i"eb.img 

echo  "Processing  EM  voxel" 

em3d  -b  -i  phaMIPG"$i".scn  -f  test  -s  2  -n  3  -t  phaLengths.AttenF.scn 

cat  testhist  »  MIPGrecon.log 

echo  ‘hostname*  »  MIPGrecon.log 

complmgs  phaMIPG"$i".img  testimg  »  MIPGrecon.log 

cp  testimg  MIPG"$i"ev.img 

end 

alias  cp  ’cp  -i’ 

File:  evai.bat 

#!/bin/csh 

#  Evaluation  of  several  reconstructed  images  using  4  reconstruction  methods 

#  File:  evai.bat  SF  20Dec93 

#  /etcA-enice  10  -u  furuie 

#  evaluate 

foreach  i  (  010  011  012  013  014  015  ) 
echo  "Processing  case""$i" 

eval3d  phaMIPG"$i".img  MIPG"$i"ab.img  testartblo.eval  4  2 
eval3d  phaMIPG"$i".img  MIPG"$i"av.img  testartvox.eval  4  2 
eval3d  phaMIPG"$i".img  MIPG"$i"eb.img  testemblo.eval  4  2 
eval3d  phaMIPG"$i".img  MIPG"$i"ev.img  testemvox.eval  4  2 
end 

File:  student.bat 

#  file:  studentbat  SF 

#  Comparison  of  reconstruction  methods  based  on  t-test 
Vm  testt 

echo  "  ****ATTENTION:  ***♦*"  >  testt 
cat  testartblo.eval  I  egrep  "Detectability"  »  testt 

#  aitblo  vs.  artvox 

echo  "♦**  artblo  vs.  artvox  ***"  »  testt 

echo  "  Table  of  t  statistics:  testartblo.eval.tbl  VS  testartvox.eval.tbl"  »  testt 

echo  "  —  HSD - CSD - SA  — "  »  testt 

echo  ""  »  testt 

echo  -n  "global  "  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  3  0  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  4  0  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  2  0  »  testt 
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echo  ""  »  testl 

echo  -n  "avCTage  "  »  testt 

student  testartblo.eval.tbl  lestartvox.eval.tbl  6  16  14  0  »  tesLt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  15  0  »  testt 
echo  ""  »  test.t 
echo  -n  "[4.0, 6.4]  "  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  8  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  11  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  5  0  »  testt 
echo  »  testt 
echo  -n  "[8.0,10.]  "  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  9  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  12  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  6  0  »  testt 
echo  ""  »  testt 
echo  -n  "[16.,20.]  "  »  testt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  10  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  13  0  »  testt 
student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  7  0  »  testt 
echo  "  "  »  testt 
echo  "  "  »  testt 

#  artblo  vs.  emblo 

echo  artblo  vs.  emblo  ***"  »  testt 

echo  "  Table  of  t  statistics:  testartblo.eval.tbl  VS  testemblo.eval.tbr  »  testt 

echo  "  —  HSD - CSD - SA  — "  »  testt 

echo  ""  »  testt 

echo  -n  "global  "  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  3  0  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  4  0  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  2  0  »  testt 

echo  ""  »  testt 
echo  -n  "average  "  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  14  0  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  15  0  »  testt 

echo  ""  »  testt 
echo  -n  "[4.0.6.4]  "  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  8  0  »  testt 
student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  11  0  »  testt 
student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  5  0  »  testt 
echo  ""  »  testt 
echo  -n  "[8.0,10.]  "  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  9  0  »  test.t 
student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  12  0  »  testt 
student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  6  0  »  tesLt 
echo  ""  »  testt 
echo  -n  "[16.,20.]  "  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  10  0  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  13  0  »  testt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  7  0  »  testt 

echo  "  "  »  testt 
echo  "  "  »  testt 

#  artblo  vs.  emvox 

echo  "***  artblo  vs.  emvox  »  testt 

echo  "  Table  of  t  statistics:  testartblo.eval.tbl  VS  testemvox.eval.tbl"  »  testt 

echo  "  —  HSD - CSD - SA  — "  »  testt 

echo  ""  »  testt 

echo  -n  "global  "  »  testt 

student  iestartblo.eval.tbl  testemvox.eval.tbl  6  16  3  0  »  testt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  4  0  »  testt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  testt 
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echo  **"  »  tesLt 

echo  -n  "average  "  »  tesLt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  14  0  »  testt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[4.0, 6.4]  "  »  tesLt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  8  0  »  tesLt 
student  testartblo.eval.tbl  teslemvox.eval.tbl  6  16  11  0  »  tesLt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  5  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[8.0,10.]  "  »  tesLt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  9  0  »  tesLt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  12  0  »  tesLt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  6  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[16„20.]  "  »  tesLt 

student  testartblo.eval.tbl  lestemvox,eval.tbl  6  16  10  0  »  tesLt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  13  0  »  tesLt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  7  0  »  tesLt 
echo  "  "  »  tesLt 
echo  "  "  »  tesLt 

#  artvox  vs.  emblo 

echo  artvox  vs.  emblo  ***"  »  tesLt 

echo  "  Table  of  t  statistics:  testartvox.eval.tbl  VS  lestemblo.eval.tbl"  »  lesLt 

echo  "  —  HSD - CSD - SA  — "  »  tesLt 

echo  ""  »  tesLt 

echo  -n  "global  "  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  3  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  4  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  2  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "average  "  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  14  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  15  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[4.0, 6.4]  "  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  8  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  11  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  5  0  »  tesLt 
echo  ""  »  lesLt 
echo  -n  "[8.0,10.]  "  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  9  0  »  tesLt 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  12  0  »  lest.t 
student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  6  0  »  test.t 
echo  ""  »  tesLt 
echo  -n  "[16.,20.]  "  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  10  0  »  tesLt 

student  testartvox.eval.tbl  lestemblo.eval.tbl  6  16  13  0  »  test.t 

student  testartvox.eval.tbl  testemblo.eval.tbl  6  16  7  0  »  tesLl 

echo  "  "  »  tesLt 
echo  "  "  »  tesLt 

#  artvox  vs.  emvox 

echo  "***  artvox  vs.  emvox  ***"  »  tesLt 

echo  "  Table  of  t  statistics:  testartvox.eval.tbl  VS  testemvox.eval.tbl"  »  tesLt 

echo  "  —  HSD - CSD - SA  — "  »  tesLt 

echo  ""  »  tesLt 

echo  -n  "global  "  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  3  0  »  tesLt 

student  testartvox.cvaI.tbl  testemvox.eval.tbl  6  16  4  0  »  lesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  tesLt 
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echo  »  test.t 

echo  -n  "average  "  »  testt 

student  testartvox.evaI.tbl  testemvox.eval.tbl  6  16  14  0  »  tesLi 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  iest.t 
echo  ""  »  tesLt 
echo  -n  "[4.0,6.4]  "  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  8  0  »  tesLt 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  11  0  »  tesLt 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  5  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[8.0,10.]  "  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  9  0  »  tesLt 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  12  0  »  tesLt 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  6  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[16., 20.]  "  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  10  0  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  13  0  »  tesLt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  7  0  »  tesLt 

echo  "  "  »  tesLt 
echo  "  "  »  tesLt 

#  emblo  vs.  emvox 

echo  "***  emblo  vs.  emvox  **♦"  »  tesLt 

echo  "  Table  of  t  statistics:  testemblo.eval.tbl  VS  testemvox.eval.tbl"  »  tesLt 

echo  "  —  HSD - CSD - SA  — "  »  tesLt 

echo  ""  »  tesLt 

echo  -n  "global  "  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  3  0  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  4  0  »  tesLl 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  tesLt 

echo  ""  »  tesLt 
echo  -n  "average  "  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  14  0  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  tesLt 

echo  ""  »  tesLt 
echo  -n  "[4.0.6.4]  "  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  8  0  »  test.t 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  11  0  »  tesLt 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  5  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[8.0,10.]  "  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  9  0  »  tesLt 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  12  0  »  tesLt 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  6  0  »  tesLt 
echo  ""  »  tesLt 
echo  -n  "[16.,20.]  "  »  tesLt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  10  0  »  tesLt 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  13  0  »  tesLt 
student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  7  0  »  test.t 
echo  "  "  »  tesLt 
echo  "  "  »  lesLt 

#  creating  table 
Nrm  tesLtbt 

echo  "  ****ATTENnON:  using  ROC  area  AVERAGE  as  FOM  for  detectability**^**"  >  tesLtbt 

echo  "  Table  output  in  file  tesLtbt" 

echo  "  artblo  artvox  emblo  emvox  "  »  tesLtbt 

#artblo 

echo  ""  »  tesLtbt 

echo  -n  "artblo  SA  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 
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student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  2  0  »  tesLtbt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  2  0  »  lesLtbt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  testtbt 

echo  »  testibt 

echo  -n  "  HSD”  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  14  0  »  tesLtbt 

student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  14  0  »  tesLtbt 

student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  14  0  »  tesLtbt 

echo  ""  »  tesLtbt 

echo  -n  "  CSD"  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

student  testartblo.eval.tbl  testartvox.eval.tbl  6  16  15  0  »  lesLtbt 
student  testartblo.eval.tbl  testemblo.eval.tbl  6  16  15  0  »  tesLtbt 
student  testartblo.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  tesLtbt 
echo  »  tesLtl^ 

#artvox 

echo  »  tesLtbt 

echo  -n  "artvox  SA  "  »  tesLtbt 

echo  -n  ”  "  »  tesLtbt 

echo  -n  **  "  >>  tesLtbt 

student  testartvox.eval.tbl  testemblo.eval.tbl  6  16  2  0  »  tesLtbt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  tesLtbt 

echo  •*"  »  tesLtbt 

echo  -n  ”  HSD"  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  -n  "  **  »  tesLtbt 

student  testartvox.eval.tbl  testemblo.eval.tbl  6  16  14  0  »  tesLtbt 

student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  14  0  »  tesLtbt 

echo  »  tesLtbt 

echo  -n  "  CSD"  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

student  testartvox.eval.tbl  testemblo.eval.tbl  6  16  15  0  »  tesLtbt 
student  testartvox.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  tesLtbt 
echo  ""  »  tesLtbt 
#emblo 

echo  ""  »  lesLtbt 

echo  -n  "emblo  SA  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  -n  ”  "  »  tesLtbt 

echo  -n  ”  "  »  tesLtbt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  2  0  »  tesLtbt 

echo  ""  »  tesLtbt 

echo  -n  "  HSD"  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  14  0  »  tesLtbt 

echo  ""  »  tesLtbt 

echo  -n  "  CSD"  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

echo  ~n  "  "  »  tesLtbt 

echo  -n  "  "  »  tesLtbt 

student  testemblo.eval.tbl  testemvox.eval.tbl  6  16  15  0  »  test.tbt 
echo  ""  »  tesLtbt 

The  next  two  tables  are  examples  of  output  of  script  student.bat. 
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Tabk  2  Table  of  t-statlstic  and  level  of  significance  associated  with  the  comparison  of  pairs  of  five  reconstruction  techniques  (example  of  output  of  script  student.bat),  using 
12  phantoms  (phaMIPG010.^01S,  phaMIPG110...115)  and  equation  (4)  for  detectability  index  and  a  summary  table  considering  HSD  and  CSD  as  average  of  ROC  area. 
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Appendix  G  FOMs  similarity 


Program  eval3d  gives  two  option  for  detectability  index  formula  (equation  (4)  and  equation  (5) 
[see  section  3.3]).  For  each  formula,  it  provides  ROC  area  (FOM,  detectability)  values  per  feature 
size  (small,  medium,  large),  average  of  these  three  ROC  areas  (_avg),  and  a  global  ROC  area 
(-SO  (see  2.4.3).  However,  there  are  some  similarities  among  these  FOMs.  Figures  10  and  11 
show  the  rank-ordering  similarity  (Yeung  and  Herman,  1989)  for  hot  spot  detectability  and  cold 
spot  detectability,  respectively.  For  a  defined  FOM,  after  a  statistical  comparison  (student  test) 
between  two  reconstruction  techniques  over  a  set  of  phantoms  (12  phantoms:  phaMIPG010..015, 
phaMIPG110...115),  we  have  a  statistic  (t).  In  our  case,  we  had  10  comparisons  (fbp,  artblob,  artvox, 
emblob,  emvox).  Now  computing  the  rank-ordering  similarity  between  two  specified  FOMs  using 
t-values,  we  have  an  entry  in  the  tables. 


Rank-ordering  Similarity  for  Hot  Spot  Detectability 

Detectability  index:  equation  5 

Detectability  index:  equation  4 

glob,  avg  small  med. 

large 

glob,  avg  small  med 

large 

0.96  0.96 

0.80  0.92 

0.68 
^ ^ 

0.68 

0.68 

0.72 

- 

0.96 

0.92 

0.92 

"I  0.9*2  ^ 

0.72 

^ - - - 

0.72 

_ 1Q0  , 

^  0^88  , 

0.96 

0.88 

^  - - — ^ 

0.68 

^  0.68 

^ - — - - 

◄ - — - 

0.88 

0.92 

- ► 

Figure  10.  Rank-ordering  similarity  for  Hot  Spot  Detectability 
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Rank-ordering  Similarity  for  Cold  Spot  Detectability 
Detectability  index:  equation  5  Detectability  index:  equation  4 

glob,  avg  small  med.  large  glob,  avg  small  med  large 


0.88  0.88  0.84  0.84 


Figure  11.  Rank-ordering  similarity  for  Cold  Spot  Detectability 

The  expected  value  and  standard  deviation  of  similarity  variable  is  (for  N=10)  0.34  and  0.1426 
respectively.  Therefore,  we  can  conclude  that  there  are  statistically  significant  similarity  between 
several  FOMs,  in  special  between  average  of  ROC  areas  (avg)  using  equation  5  and  4;  average  of 
ROC  areas  and  ROC  area  of  small  features;  and  global  ROC  area  and  ROC  area  of  small  features.  The 
tables  also  show  that  small  features  have  stronger  influence  on  the  calculation  of  global  and  average 
ROC  areas  than  medium  and  large  sizes  for  the  studied  phantoms. 
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